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Single  crystal  silicon  was  selected  as  a  model  material  in  which  to  study  the 
correlation  of  the  fracture  surface  features  as  characterized  by  their  fractal  dimension  for 
different  orientations  of  fracture  with  the  fracture  toughness  of  the  material  as  measured 
using  the  strength-indentation,  and  fracture  surface  analysis  techniques.  Single  crystal  Si 
was  selected  for  several  reasons:  Si  is  a  brittle,  monoatomic  material  that  will  obviate  the 
complication  of  microstructure,  i.e.,  grains,  pores,  etc.  And  Si  has  been  well  studied  so 
that  many  of  its  properties  are  well  characterized.  Flexure  bars  were  indented  with  a 
Vickers  indent  at  various  loads  and  fractured.  After  calculating  the  fracture  strength  and 
toughness,  the  surfaces  were  analyzed  and  characterized  using  fracture  surface  analysis 
and  slit-island  analysis.  These  analyses  provided  the  size  of  the  fracture  initiating  defect, 
the  geometry  of  the  surrounding  topography  including  the  location  of  the  regions  of  crack 
branching  and  the  fractal  dimensions  of  selected  areas  on  the  fracture  surface. 

The  {100}  plane  is  found  experimentally  to  be  the  one  with  the  highest  toughness 
of  the  three  studied  planes.  The  fracture  plane  on  the  {100}  and  {110}  fracture  plane 
both  initiate  on  the  original  plane  and  have  the  tendency  to  deviate  to  the  {111}  plane. 


The  fracture  surfaces  of  Si  have  been  analyzed  and  are  found  to  be  fractal  both  at 
the  atomic  scale  using  the  scanning  tunneling  microscope  in  a  previous  study  and  at  the 
micrometer  scale  using  the  optical  microscope  in  this  study.  The  irregularity  of  the 
fracture  surface  is  too  complicated  for  Euclidean  geometry  to  describe  easily.  Self- 
similarity  and  scale  invariance  is  suggested  by  the  fracture  surface  appearance.  The 
fractal  dimension  is  higher  for  the  fracture  plane  with  the  higher  toughness  value.  It  is  in 
agreement  with  other  studies  which  found  higher  toughness  values  associated  with  higher 
fractal  dimensions. 

We  also  demonstrate  the  formation  of  fractal  surfaces  during  fracture  using  a 
molecular  dynamics  (MD)  approach  in  a  single  crystal  silicon  structure  and  compare  the 
simulated  results  with  experimental  work.  MD  simulations  using  Stillinger- Weber 
potential  and  Coulombic  (Modified  Bom-Mayer)  potential  are  performed  in  different 
orientations  to  investigate  the  difference  in  toughness  and  fractal  dimension  with  respect 
to  orientations.  The  close  agreement  between  simulated  and  measured  fractal  structures 
of  the  fracture  surface  suggests  that  this  is  a  promising  method  for  investigating  atomic- 
level  processes  during  fracture. 
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CHAPTER  1 
INTRODUCTION 


The  strength  of  brittle  materials  is  determined  by  the  stress  required  to 
spontaneously  propagate  a  crack.  The  strength  is  dependent  on  the  local  stress 
distribution  around  a  crack  tip  which  is  often  treated  using  a  stress-intensity  factor 
approach.  The  fracture  toughness,  which  is  a  measure  of  the  resistance  to  crack 
growth,  defines  the  work  or  mechanical  energy  expended  in  propagating  the  crack. 
The  fracture  toughness  can  be  represented  by  the  critical  stress  intensity  factor,  the 
fracture  energy  or  the  critical  strain  energy  release  rate.  Techniques  to  determine  the 
fracture  toughness  of  ceramic  materials  include  cleavage,  double  cantilever  beam 
(DCB),^  crack  indentation,''^  strength-indentation,^  and  fracture  surface  analysis  (FSA)^ 
techniques.  The  mechanism  of  fracture  in  brittle  materials  is  strongly  influenced  by 
processes  at  the  atomic  level.  The  subject  of  fracture  in  brittle  materials  has  long  been 
of  interest;  however,  the  underlying  processes  have  remained  difficult  to  investigate 
directly.  Some  factors  such  as  phonon-phonon  interactions  have  been  analyzed 
recently,  but  the  dynamics  of  bond  breaking  in  brittle  materials  is  only  beginning  to  be 
studied. ^"^^  Although  it  is  intuitively  obvious  that  there  must  be  a  relationship  between 
the  bond  breaking  process  and  the  macroscopic  measures  of  fracture  such  as  fracture 
toughness  and  fracture  topography,  there  has  been  little  quantitative  evidence  of  the 
connection.  Fractal  geometry,  which  is  a  non-Euclidean  geometry  exhibiting  self- 
similarity  and  scale  invariance,  is  a  new  tool  that  can  be  used  to  relate  atomic  scale 
processes  to  macroscopic  processes  and  features.   The  fracture  surfaces  of  Si  have  been 


analyzed  and  are  found  to  be  fractal  both  at  the  atomic  scale  using  the  scanning 
tunneling  microscope  and  at  the  micrometer  scale ^^  using  the  optical  microscope. 

Fracture  markings  on  glasses  and  polycrystalline  brittle  materials,  known  as 
mirror,  mist,  and  hackle,  are  precursors  to  crack  branching  and  can  be  used  to 
describe  the  stress  state^  and  characteristics  of  crack  propagation.^^  These  marking 
have  been  observed  for  more  than  50  years  and  were  related  quantitatively  to  the  stress 
condition  in  the  1950's.^^  More  recently,  the  repetition  of  these  features  was  observed 
and  quantitatively  related  to  stress  intensity.  Ravi-Chandar  and  Knauss^^  also  noted 
that  mist  and  hackle  are  self-similar;  i.e.,  they  appear  to  be  physically  similar  and 
produced  in  the  same  fashion.  Their  description,  however,  did  not  emphasize  the  self- 
similar  namre  of  the  features.  The  distance  to  branching  is  directly  related  to  the  size 
of  the  fracture  initiating  crack.  Thus,  one  can  imagine  that  the  fractore  surface  features 
can  be  related  to  the  bond  breaking  processes  at  fractures. 

Mandelbrot^  ^  synthesized  a  branch  of  mathematics  that  provides  a  tool  with 
which  to  analyze  self-similar  and  self-affine  processes.  A  self-similar  process  is  one  in 
which  a  feature  at  one  magnification  is  related  to  another  at  another  magnification  by  a 
scalar  quantity.  A  self-affine  process  is  one  in  which  this  magnification  factor  is  a 
vector  quantity.   Mandelbrot  called  this  branch  of  mathematics  "fractal  geometry. " 

Mecholsky  et  al.^^  summarized  that  fracture  in  brittle  materials  can  be  described 
as  a  self-similar  process  and  thus,  can  be  mathematically  described  by  fractal  geometry. 
They  derived  an  expression  which  relates  the  fracmre  of  bonds  at  the  atomic  level  (y)  to 
a  characterization  of  the  resulting  morphology  on  the  fracture  surface  of  polycrystalline 
materials,  D*,  which  is  the  fractional  part  of  the  fractal  dimension.  They  also  derived 
a  relationship  between  the  surface  energy  and  the  fracture  energy  of  single  crystals, 
polycrystalline  ceramics,  and  glass-ceramics.  Their  results  suggest  that  modeling 
brittle  fracture  as  a  fractal  process  may  be  useful  in  distinguishing  toughening 
mechanisms  and  in  relating  atomic  bonding  and  fracture-surface  morphology.     The 


present  research  was  undertaken  to  study  the  relationship  between  fracture  surface 
morphology  of  single  crystals  and  their  resistance  to  fracture. 

Unlike  mathematically  generated  fractals,  the  self-similarity  of  fracture  surface 
geometry  is  bounded  by  measurement  limitations.  These  limitations  will  be  discussed 
later.  However,  the  self-similar  nature  of  fractal  surfaces  offers  a  means  of  scaling 
macro-scopically  analyzable  structures  to  microscopic  processes  on  the  atomic  scale. 
In  addition,   the  fractal  dimension  has  been  shown  to  relate  dkectly  to  fracture 

toughness,  thus  forming  a  link  between  average  macroscopic  behavior  and  possible 
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atomic  processes.    ' 

Single  crystal  silicon  was  selected  as  a  model  material  in  which  to  study  the 
correlation  of  the  fracture  surface  features  as  characterized  by  their  fractal  dimension 
for  different  orientations  of  fracture  with  the  fracture  toughness  of  the  material  as 
measured  using  the  strength-indentation,  and  fracture  surface  analysis  techniques. 
Single  crystal  Si  was  selected  for  several  reasons:  Si  is  a  brittle,  monoatomic  material 
that  will  obviate  the  complication  of  microstructure,  i.e.  grains,  pores,  etc.  And  Si  has 
been  well  studied  so  that  many  of  its  properties  are  well  characterized.  Flexure  bars 
were  indented  with  a  Vickers  indent  at  various  loads  and  fractured.  After  calculating 
the  fracture  strength  and  toughness,  the  surfaces  were  analyzed  and  characterized  using 
fracture  surface  analysis  and  slit-island  analysis.  These  analyses  provided  the  size  of 
the  fracture  initiating  defect,  the  geometry  of  the  surrounding  topography  including  the 
location  of  the  regions  of  crack  branching  and  the  fractal  dimensions  of  selected  areas 
on  the  fracture  surface. 

In  this  study,  we  demonstrate  the  formation  of  fractal  surfaces  during  fracture 
using  a  molecular  dynamics  (MD)  approach  in  a  single  crystal  silicon  structure  and  try 
to  compare  the  simulated  results  with  the  experimental  works.  The  close  agreement 
between  simulated  and  measured  fractal  structures  of  the  fracture  surface  suggests  that 
this  is  a  promising  method  for  investigating  atomic-level  processes  during  fracture. 


The  use  of  molecular  dynamics  for  studying  the  failure  of  materials  has  recently 
become  of  great  interest  to  investigators.  Works  by  Paskin  et  al.,l9  Dienes  and 
Paskin,20  and  other  collaborative  works  involving  Paskin,  have  performed  simulations 
of  crack  growth  in  two-dimensional  lattices.  These  works  demonstrated  the 
applicability  of  the  Griffith  criteria  in  simple  2-D  systems  modeled  by  molecular 
dynamics,  and  studied  mechanisms  of  crack  growth,  including  the  formation  of 
dislocations  and  speed  of  crack  propagation,  as  functions  of  applied  stress,  crack 
length,  and  potential  parameters.  These  works  provide  an  atomistic  basis  to  fracture 
mechanics;  however,  they  do  not  provide  details  on  the  fracture  of  particular  materials. 

Soules  and  Busby  used  molecular  dynamics  to  study  the  rheological  properties 
and  fracture  behavior  of  a  sodium  silicate  glass.  Among  die  experiments  in  this  work, 
a  uniaxial  tension  was  applied  to  a  system  and  the  resulting  behavior  was  studied.  The 
periodic  boundary  conditions  used,  i.e.,  on  only  two  sides,  made  this  experiment 
equivalent  to  the  fracture  of  a  50A  diameter  fiber,  so  bulk  properties  were  not  studied, 
nor  even  approximated.  This  work  was  followed  by  works  of  Soules,^  Ochoa  and 
Simmons,  and  Ochoa  et  al.,  which  demonstrated  the  importance  of  dynamic  effects 
in  the  fracture  of  solids.  Works  by  Ochoa,  et  al.,^*^  made  use  of  full  periodic  boundary 
conditions,  i.e.  on  six  sides,  thus  approximating  fracture  within  a  bulk  material.  These 
simulations  collectively  studied  the  overall  response  of  a  system  to  an  applied  strain, 
but  did  not  study  the  atom-by-atom  process  in  detail.  A  work  by  Simmons  et  al.^^ 
examined  the  mechanisms  involved  in  an  individual  bond  fracture  to  describe  a  process 
by  which  an  entire  fracture  surface  may  be  stabilized.  Thus,  these  studies  give  possible 
details  on  the  fracture  process  for  silicate  glasses,  and  infer  details  of  fracture  for  other 
amorphous  solids. 

The  difference  in  fracture  between  amorphous  silica  and  single  crystal  silicon  is 
that  silica  has  no  preferred  cleavage  plane  while  single  crystal  silicon  prefers  fracturing 
on  the  {111}  and  {110}  plane. -^^    Thus,  different  loading  orientations  under  different 


loading  histories  in  single  crystal  silicon  would  result  in  different  fracture 
topographies^'^''^^  (although  the  fracture  processes  may  be  the  same  for  different 
orientations).  The  effect  of  elastic  anisotropy  on  fracture  provided  the  inspiration  to 
use  single  crystal  silicon  for  the  molecular  dynamics  simulation.  From  the  simulated 
loading  history,  the  comparison  of  fracmre  toughness  in  different  orientations  can  be 
made.  From  the  simulated  fracture  topography,  the  irregularity  of  the  fracture  surface 
in  different  orientations  can  be  analyzed  using  fractal  analysis  and  the  fractal  dimension 
can  be  compared  with  the  toughness.  Kieffer  and  Angell^^  used  the  MD  simulation 
method  to  generate  stable  silica  aggregates  at  various  low  densities  sunilar  to  those  of 
experimental  aerogels.  They  found  that  fractal  dimensions  and  range  of  self-similarity 
can  be  extracted  from  the  radial  distribution  functions  ui  those  structure  obtained  from 
MD  simulation.  They  believed  that  a  non-integer  dimension  is  a  characteristic  feature 
of  the  aerogel  structure  used  in  their  smdy.  However,  because  the  way  which  they 
used  the  change  of  the  slope  of  the  radial  distribution  functions  with  respect  to  the 
density  to  obtain  fractal  dimension  in  their  study  was  not  reasonable,  thus  they  got  a 
weird  result  which  shows  that  a  denser  aerogel  has  a  lower  fractal  dimension.  It  is 
contradict  to  the  fact  that  in  the  nature  a  full  occupied  volume  has  a  dimension  3  while 
a  volume  with  porosity  has  dimension  less  than  3.  Thus  we  must  use  a  different 
approach  from  theirs  to  obtained  the  fractal  dimension  in  this  study.  Here  the  fractal 
dimension  on  the  generated  fracture  surfaces  at  the  atomic  scale  from  MD  simulation 
will  be  analyzed  using  slit-island  technique  and  Richardson-plot. 

This  investigation  will  address  several  topics.  First,  there  will  be  an 
experimental  determination  of  the  fracmre  toughness  in  single  crystal  silicon  as  a 
function  of  orientation  of  the  crystal  plane.  Second,  the  fractal  dimension  of  the 
fracmre  surfaces  of  the  single  crystal  silicon  for  the  different  orientations  will  be 
determined  experimentally.  The  fractal  dimensions  will  be  determined  using 
Richardson  plots  and  slit-island  analyses  of  contours  of  (horizontal)  sections  through 


the  fracture  surfaces.  Third,  MD  simulations  using  Stillinger- Weber  and  Coulombic 
potentials  will  be  performed  in  different  orientations  to  investigate  the  difference  in 
toughness  and  fractal  dimension  with  respect  to  orientations.  Finally,  there  will  be 
discussions  about  the  relationship  between  the  fracture  toughness  and  the  fractal 
dimensions  of  single  crystal  silicon  for  the  different  orientations  and  also  comparisons 
between  experimental  work  and  simulation  results. 


CHAPTER  2 
FUNDAMENTALS 


Structure  of  Single  Crystal  Silicon 

Silicon  is  bonded  by  four  covalent  bonds  which  produce  a  tetrahedron.  Each 
silicon  atom  is  bonded  with  four  other  silicon  atoms  due  to  the  nature  of  the  covalent 
bonding.  The  covalent  bond  formed  by  two  Si  atoms  sharing  electrons  is  very 
localized  and  directional.  The  space  lattice  of  the  diamond  structure  is  face-centered 
cubic  with  two  atoms  per  lattice  site,  one  at  0,0,0  and  the  other  at  1/4,  1/4,  1/4  (Fig. 
2.1).  Each  atom  is  tetrahedrally  bonded  to  four  nearest-neighbors  due  to  the  nature  of 
covalent  bonding.  As  these  tetrahedral  groups  are  combined,  a  large  cube  is 
constructed,  which  is  called  a  diamond  cubic  (DC)  unit  cell.  The  large  cube  contains 
eight  smaller  cubes  that  are  the  size  of  the  tetrahedral  cube  but  only  four  of  the  cubes 
contain  tetrahedrons.  The  lattice  is  a  special  face-centered  cubic  (FCC)  structure, 
which  is  shown  in  Fig.  2.1.  The  atoms  on  the  corners  of  the  tetrahedral  cubes  provide 
atoms  at  each  of  the  regular  FCC  lattice  points.  Four  additional  atoms  are  present 
within  the  DC  unit  cell  from  the  atoms  in  the  center  of  the  tetrahedral  cubes. 
Therefore,  there  are  eight  atoms  per  unit  cell.  The  lattice  parameter  length  is  5.43  A 
and  the  unit  cell  diagonal  length  is  9.41  A  as  shown  in  Fig.  2.1.  The  Si  bond  distance 
is  2.34  A. 

Single  crystal  Si  is  generally  classified  as  a  brittle  material  whose  atoms  are 
strongly  covalently  bonded.  The  primary  cleavage  plane  is  the  {111}  plane.  The 
principal  factor  for  crack  initiation  is  bond  rupture. 


(a) 
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(c) 


Fig.  2.1  (a)  Tetrahedron  and  (b)  the  diamond  cubic  (DC)  unit  cell,  (c)  The  length 

of  AB  is  5.43  A,  AC  is  9.41  A,  and  AD  is  7.68  A. 
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Fracture  of  single  crystal  Si  on  the  {100},  {110}  and  {111}  plane  has  been 
observed. ^^'^^  Not  only  the  toughness,  but  also  the  fracture  surface  topography  were 
found  to  be  different  on  the  different  planes.  These  observations  mean  that  elastic 
anisotropy  plays  an  important  role  in  fracture;  the  bond  breaking  process  is  considered 
the  same  since  only  one  kind  of  Si-Si  covalent  bond  exists. 

Failure  Analysis 

The  characterization  of  flaws  and  their  relationship  to  strength  plays  an 
important  roll  in  the  analysis  of  brittle  materials.  By  comparing  fracture  mechanics 
relationships  with  observed  fracture  surface  markings,  important  information  about  the 
fracmre  toughness  can  be  determined.  This  section  will  present  the  basic  principles  of 
fracture  mechanics  and  fracture  surface  analysis  of  brittle  materials. 

Fracture  Mechanics 

One  of  the  first  analysis  of  fracture  behavior  of  components  that  contain  sharp 
discontinuities  was  developed  by  Griffith.^^'^^  This  analysis  was  based  on  the 
assumption  that  incipient  fracture  in  ideally  brittle  materials  occurs  when  the  magnimde 
of  the  elastic  energy  supplied  at  the  crack  tip  during  an  incremental  increase  in  crack 
length  is  equal  to  or  greater  than  the  magnitude  of  the  elastic  energy  at  the  crack  tip 
during  an  incremental  increase  in  crack  length.  By  using  a  stress  analysis  developed  by 
Inglis,  Griffith  related  fracture  stress  to  the  flaw  size  by  means  of  an  energy  balance. 
The  relationship  can  be  written  as: 


rov'.,\^'^ 


f^f 


2Ey 


(2.1) 


where  af  is  the  failure  stress. 
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E  is  the  modulus  of  elasticity, 

Y  is  specific  surface  energy,  and 

c  is  the  critical  flaw  size. 
Although  Griffith  was  the  first  to  analyze  the  relationship  between  strength  and 
flaw  size,  it  was  Irwin  who  developed  fracture  mechanics  into  the  present  day  form. 
Irwin  followed  the  work  done  by  Griffith,  Orowan,  and  Inglis  to  develop  what 
is  known  as  linear  elastic  fracture  mechanics.  Irwin  first  analyzed  the  fracture  of 
flawed  components  using  a  stress  analysis  based  on  the  Westergaard^^  solution  of  an 
elliptical  crack  in  an  infinite  plate.  For  a  surface  cracked  specimen  under  mode  I 
(tensile)  loading,  Irwin  derived  the  following: 

CTf  =  -^  ,  (2.2) 


where  of  is  the  failure  stress, 

Kjc  is  the  critical  stress  intensity  factor  for  mode-I  loading, 

B  is  a  geometrical  factor  which  accounts  for  flaw  shape,  location, 

and  loading  geometry,  and 

c  is  the  critical  flaw  size. 
Irwin  also  analyzed  fracture  in  terms  of  the  strain  energy  release  rate,  G.    G  is 
defined  as  the  elastic  energy  per  unit  crack  length,  U/c,  and  can  be  related  to  the 
failure  stress: 


af  = 


^EG.V^^ 


'c 


(2.3) 


y  nc  J 
where  G^  is  the  elastic  energy  release  rate  at  fracture. 

The  strain  energy  release  rate  can  be  related  to  surface  energy,  y,  and  it  is  as 
follows  for  plane-stress  condition: 

G  =  2y.  (2.4) 
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In  the  above  equation,  only  surface  energy,  y,  is  referred.  Indeed  during  fracture 
process,  the  released  strain  energy  not  only  produce  two  free  surface,  but  also  produce 
different  terms  like  sound,  light,  electron  emission.  Thus,  the  y  used  here  indeed  is  not 
just  as  simple  as  surface  energy  and  some  other  energy  is  also  included. 

The  stress  intensity  factor  can  be  related  to  the  strain  energy  release  rate  and  the 
fracture  surface  energy.  By  combining  those  previous  equations,  it  can  be  shown  that 


Kic=VEG=V2yE.  (2-5) 

Fracture  Surface  Analysis 

Useful  information  on  the  mechanical  and  fracture  behavior  of  a  failed  ceramic 
can  be  gained  from  microscopic  examination  of  the  fracture  surface.  "  The  fracture 
surface  of  brittle  materials  exhibits  distinct  fracture  markings  surrounding  the  critical 
flaw.  These  features  shown  in  Figure  2.2  can  be  related  by  known  fracture  mechanics 
relationships  to  provide  additional  information  on  the  fracture  process  and  can  also  be 
used  to  locate  the  origin  of  failure. 

The  fracture  surface  markings  shown  m  Figure  2.2  can  be  used  to  describe  the 
stress  state  of  a  brittle  material.  ^"^'^^'-^^'^^  Four  different  regions  referred  to  as  the 
fracture  mirror,  mist,  hackle  and  crack  branching  can  be  seen  in  Figure  2.2.  These 
regions  are  associated  with  particular  stress  intensity  levels  and  crack  velocities  which 
are  responsible  for  each  distinct  region.^^  As  a  crack  propagates  from  the  critical  flaw, 
a  smooth  region  which  is  basically  perpendicular  to  the  tensile  axis  is  formed.  This 
smooth  region  is  known  as  the  fracture  mirror.  The  fracture  mirror  is  typically 
bounded  by  a  region  of  small  radial  ridges  known  as  mist  and  the  mist  is  bounded  by 
an  even  rougher  radial  ridge  region  called  hackle.  Finally,  the  propagating  crack 
reaches  a  characteristic  energy  level  and  crack  branching  occurs.  Crack  branching  is  a 
region  where  two  or  more  cracks  form  from  the  primary  crack  front. 
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»    k4Sl^-J 


■4^. 


Fig.  2.2  Schematic  of  features  found  on  the  fracture  surface  of  a  brittle  material 

subjected  to  a  constant  load.  Solid  semi-elliptical  line  at  center  represents  the  initial 
flaw  size  depth  aj  and  width  2bi.  Dashed  line  represents  the  outline  of  critical  flaw 
depth  a^r  and  width  2bci-.  Mirror/mist,  r^,  mist/hackle,  T2,  and  crack  branching,  x^, 
radii  are  shown  along  the  tensile  axis. 
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Once  the  failure  origin  is  located  after  fracture,  its  size  can  be  measured  and 
used  to  determined  the  toughness  of  the  material.  Figure  2.3  shows  a  semi-elliptical 
crack  in  a  plate.    It  can  be  shown  through  a  Griff ith/Irwin  approach  that  the  stress 

intensity  on  a  semi-elliptical  surface  crack  of  depth,  a,  and  half-width,  b,  can  be  given 

35  37 
as:    ' 

^^       V'l.27TaVa  .^  ,, 

Ki^ '  (2.6) 

where  a  is  the  semi-minor  axis  for  an  elliptical  crack, 

a  is  the  applied  stress,  and 

O  is  an  elliptical  integral  of  the  second  kind.  O  varies  between  <I>=  1 
for  a  slit  crack  (a/b=0)  to  O  =  1.57  for  a  semi-circular  crack  (a/b  =  l). 
The  criterion  for  failure  in  a  brittle  material  is  that  Kj  >  Kj^,  at  which  point  the 
initial  flaw  will  propagate  spontaneously,  where  Kj,  the  stress  intensity,  is  a  measure  of 
the  magnification  of  the  external  loading  at  the  crack  tip,  and  Kj^  is  the  critical  stress 
intensity  factor  or  the  fracture  toughness  of  the  material.  The  strength  or  failure  stress, 
Of,  of  a  brittle  material  can  be  related  to  the  flaw  size: 

af  =  ^£l^.  (2.7) 

•\/1.27lC 

For  a  semi-circular  surface  crack,  which  is  stress  free  and  which  is  small  relative  to  the 
thickness  of  the  cross-section,  O  =  1.57  and  Kj^  can  be  found  as: 

Kic  =  1.24afA/^,  (2.8) 

where  c  is  the  square  root  of  ab  and  is  the  size  of  a  semi-elliptical  crack  which  has  the 

30 

size  of  an  equivalent  semi-circular  crack. 

Thus,  the  fracture  toughness  of  a  material  can  be  calculated  if  the  failure  stress 
is  determined  and  the  crack  dimensions  a  and  b  are  measured. 
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Fig.  2.3  An  semi- infinite  plate  under  uniform  tension  and  containing  a  semi- 

elliptical  crack. 
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Indentation  Fracture  Mechanics 

Indentation  fracture  mechanics  has  been  well  established  as  an  important  method 
for  the  study  of  the  mechanical  behavior  of  brittle  materials. ^^"'^'^  Indentation  provides 
a  means  of  introducing  flaws  into  a  material  with  a  controlled  size,  shape,  and  location 
on  the  sample.  This  is  in  contrast  to  materials  with  namrally  occurring  flaws  where  the 
flaws  vary  in  size,  shape,  location  and  concentration.  Indentation  provides  a  means  of 
characterizing  the  fracture  process  by  introducing  a  controlled  flaw  into  a  brittle 
material.  The  following  discussion  will  concentrate  on  fracture  surface  analysis  and 
strength  indentation  analysis. 

Fractography-Indentation  Analysis 

Figure  2.4(a)  shows  a  schematic  of  the  indentation  deformation/fracture  pattern 
for  the  Vickers  diamond  geometry:  P  is  the  peak  load  and  a  and  b  are  characteristic 
dimensions  of  the  inelastic  impression  of  the  radial/median  crack,  respectively.  Figure 
2.4(b)  shows  the  schematic  of  a  fiiUy  developed  Vickers  indent.  The  flaw  is  taken  to 
possess  a  penny-like  geometry.  The  following  discusses  the  cases  of  zero  and  non  zero 
residual  stress  terms  separately: 

(1)  Zero  residual  contact  stresses:  If  the  indentation  flaw  were  to  be  free  of 
residual  stresses,  the  stress  intensity  factor  for  uniform  tensile  loading  would  have  the 
standard  form 

Ka=Yafcl/^  (2.9) 

where  K^  is  the  applied  stress  intensity  factor  without  residual  stress, 

Y  is  a  unitless  geometrical  factor.  For  a  semi-circular  crack,  Y  =  1.24, 
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Fig.  2.4  Configuration  of  median/radial  cracks  for  Vickers  indentation  showing: 

(a  normal  indent  load,  P,  generating  median  opening  forces  Pg  (elastic  component)  and 
Pj.  (residual  component),  (b)  load  removal  eliminates  the  elastic  component,  and  (c) 
fully  developed  radial  crack  pattern. 
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Qf  is  the  failure  strength,  and 
c  is  the  critical  flaw  size. 
(2)  Nonzero  residual  contact  stresses:   With  residual  contact  stresses  present  at 
the  flaw  origin  it  becomes  necessary  to  include  an  additional  tensile  term  in  the  stress 
intensity  factor: 44 

Ki  =  Ka  +  Kr  =  YaaC^^2  +  x^p  /  ^^'^ ,  (2. 10) 

where  Kj.  is  the  stress  intensity  factor  with  residual  stress, 

P  is  the  peak  load,  and 
Xj  is  a  unitless  geometry  constant. 
For  large  c  values,  the  applied-stress  term  controls  the  fracture,  as  before,  for 
small  c  values,  it  is  the  residual-stress  term  that  dominates.  "Under  equilibrium 
fracture  conditions  the  flaw  will  accordingly  undergo  a  precursor  stage  of  stable  growth 
as  the  tensile  loading  is  applied;  failure  then  occurs  when  the  crack  reaches  a  critical 
size,  at  which  point  the  applied  stress  is  intense  enough  to  cause  spontaneous 
propagation."^^  This  critical  configuration  is  obtained  by  inserting  Kj  =  K^  (Kj.  is 
used  here  instead  of  Kj,,  because,  strictly  speaking,  the  toughness  determined  using 
indentation  is  different  than  Kj^;  however,  in  practical  terms,  they  are  equivalent)  into 
the  above  equation  and  evaluating  the  instability  condition  da^/dc  =  0: 

^^r4x,pf^^r3K,  f 

'^     I   Kc    J  l4Y(7f  J  ' 

from  which  Kp  can  be  found  as: 45, 46 


(2.11) 


Ke  =  -YcjfC^^^  (2.12) 

For  a  semi-circular  flaw,  which  is  small  relative  to  the  thickness  of  the  material  into 
which  it  is  placed,  Y  =  1.24,  and  the  above  equation  can  be  rewritten  as: 

Kc  =  1.65afC^/^  (2.13) 
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Thus,  the  resistance  to  fracture  can  be  determined  from  a  measurement  of  the  surface 
crack  size  and  the  strength  of  bars  that  have  been  indented. 


Strength-Indentation  Technique^ 


Consider  the  Vickers  diamond- induced  radial  cracks,  which  are  shown  in  Fig 
2.4.  If  the  applied  loading  is  uniaxial,  the  indentation  is  aligned  with  one  set  of 
pyramidal  edges  parallel  to  the  tensile  axis. 

It  has  been  shown  that  the  toughness  may  be  derived  using  the  following 
expression: 

K.=.^(|)"^(.,P-r,  (2,4) 

,  R      I  256^^2/3f-R  I 

where  r]^  =   Y     Sy       is  a  geometry  constant, 

\  z  /  J 

E  is  the  Young's  modulus,  and 

H  is  the  hardness. 
Replacing  iq^  by  an  average  quantity,  0.59,  would  add  no  more  than  10%  to  the  error 
in  the  Kc  evaluation  for  a  material   whose  elastic/plastic  parameters   are  totally 
unknown.   So  the  above  equation  can  be  rewritten  as 


Kc=0.59 


1/8  ,2/4 

(afPl/3)      .  (2.15) 


f-1 

Thus,  the  resistance  to  fracture  can  also  be  measured  from  a  flexural  test  in  which  the 
specimens  have  been  indented  at  different  load  levels  without  measuring  the  crack  size. 
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Fractal  Geometry 

Fractal  geometry  is  a  non-Euclidean  geometry  that  was  rediscovered, 
popularized  and  applied  by  B.  B.  Mandelbrot.  ^^  The  word  fractal  is  derived  from  the 
Latin  "fractus"  which  means  fragmented  or  broken.  Nature  exhibits  diverse  structures 
having  inherent  irregularities.  Some  of  these  structures  can  be  described  using 
Euclidean  geometry  but  other  structures  are  better  described  by  using  fractal  geometry. 
In  his  book,  The  Fractal  Geometry  of  Natare.^^  Mandelbrot  describes  the  extensive 
applications  of  fractal  geometry.  The  concepts  have  been  used  to  described  the 
geometry  of  clouds,  soot  aggregates,  dielectric  breakdown,  coastlines,  and  other 
natural  phenomena,  including  fracture. 

Fractal  geometry  exhibits  self-similarity  or  self-affinity,  shows  scale  invariance, 
and  is  characterized  by  the  fractal  dimension.  Self-similarity  means  that  features  in 
different  regions  appear  to  be  similar  to  one  another.  They  can  be  related  by  a  scalar 
multiple.  Features  in  a  self-affine  object  are  related  by  a  vector  quantity.  Scale 
invariance  occurs  if  a  feature  on  different  scales  appears  to  be  the  same.  These 
properties  indicate  that  a  fractal  can  be  created  by  using  a  shape  to  be  repeated  and 
reduced  in  size  following  a  prescribed  sequence. 

Fractal  Dimension 

In  Euclidean  geometry,  objects  occupy  integer  dimensions,  i.e.,  1,  2,  3,  and  so 
on.  Fractal  geometry  admits  to  the  description  of  objects  that  occupy  fractional 
dunensions.  Fractal  geometry  is  characterized  by  its  fractal  dimension  which  is  related 
to  measurement  theory. 
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Consider  how  we  measure  the  length  of  a  line.  The  measurement  of  a  line 
requires  a  measuring  tool  or  ruler.  A  line  is  given  in  Fig.  2.5(a).  The  length  of  this 
line  is  measured  by  covering  it  completely  with  as  few  of  a  given  radius,  R,  discs  as 
possible,  as  in  Fig.  2.5(b).   The  length  is  given  by 

length  —  (number  of  covering  discs)  x  (  2R).  (2.16) 

If  we  have  a  new  supply  of  discs  with  a  different  radius,  say  R/2,  then  the  line  length  is 
determined  in  the  same  manner,  and  a  new  line  length  is  computed,  as  shown  in  Fig. 
2.5(c).  Note  the  difference  in  the  measured  line  length  in  Fig.  2.5(c).  Although  the 
measuring  disc  was  decreased  to  half  of  the  original  measure,  the  number  required  to 
cover  the  line  more  than  doubled.  The  "tortuous"  nature  of  this  curve  gives  a  length 
that  is  scale  dependent.  As  the  measuring  scale  decreases,  the  measured  length 
increases. 

The  measured  length  of  a  line  can  be  described  by  the  following  equation, 
called  Richardson's  equation: 

L  =  kS^"°,  (2.17) 

where  L  is  the  measured  length, 

S  is  the  measured  scale, 
D  is  the  dimension  of  the  curve, 
k  is  a  proportionality  constant. 

If  D  is  equal  to  one,  the  line  is  Euclidean  and  its  measured  length  is  not  a 
function  of  scale.  If  1<  D  <  2,  the  curve  is  said  to  be  fractal  with  a  dimension  given  by 
D.  Notice  that  for  a  given  line  segment,  the  length  is  not  constant,  but  the  dimension 
is.  Thus,  the  fractal  dimension  is  the  dimension  in  which  a  measure  can  be  made,  not 
the  measure  itself. 

The  fractal  dknension  of  a  line  can  be  computed  from  the  above  equation.  The 
line's  length  is  computed  over  a  range  of  scales.  A  log-log  plot  of  length  vs.  scale 
gives  a  straight  line  with  slope  equal  to  1-D,  as  shown  in  Fig.  2.6. 
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(a) 


Fig.  2.5  Measurement  of  a  line  segment  (  or  signal),   (a)  A  schematic  for  a  line 

which  has  texture.  Measuring  the  length  of  a  line  requires  a  ruler.  The  (circle) 
covering  has  a  diameter,  2R.  (b)  The  line  is  completely  covered  to  determine  its 
length,  (c)  Changing  the  size  of  the  measuring  disc  may  produce  a  different  measure  of 
the  length. 
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(b) 


O 


diameter  =  2R 


length  =  #  of  disc  x  diameter  2R 


(c) 


o 


diameter  =  2r 


length  =  #  of  disc  x  diameter  2r 


Fig.  2.5  Continued. 
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Slope  =1-D 


logL 


LogE 


Fig.  2.6  A  log-log  plot  of  length  vs.  scale  gives  a  straight  line  with  slope  equal  to 

1-D. 
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Mathematical  fractals  will  obey  the  above  equation  for  all  scales.  However, 
physical  fractals  may  have  a  finite  cutoff  at  both  large  and  small  scales.  In  the  physical 
world,  the  fluctuation  of  a  line  will  fmd  its  limit  in  the  smallest  measurable  feature. 
The  ruler  length  determines  the  smallest  feature  we  can  measure.  If  variations  in 
strucmre  are  below  the  sensitivity  of  the  measurement  scale,  the  line  will  appear  to  be 
Euclidean.  In  that  case,  the  longest  line  we  could  measure  would  depend  on  the  finest 
possible  scale  of  measurement.  We  could  not  conclude  that  the  line  has  no 
fluctuations,  only  that  we  have  no  ability  to  measure  the  flucUiations. 

On  the  contrary,  a  line  may  be  composed  of  indivisible  components.  After  our 
measuring  stick  has  become  smaller  than  the  smallest  of  these  components,  the  line 
would  cease  to  exhibit  scale-dependent  length.  This  length,  then,  would  represent  a 
true  upper  bound.  Such  a  line  could  be  called  "fractal"  over  the  range,  it  exhibited 
scale-dependent  length.  It  is  not,  however,  a  fractal  in  the  strict  mathematical  sense. 

The  combination  of  the  two  lunits  would  result  in  the  Richardson  plot  as  shown 
in  Fig.  2.7.  The  object  would  be  described  as  a  fractal  over  a  bounded  range  of  scales. 
The  cut-off  scales  are  of  special  significance  to  the  particular  geometry  of  a  specific 
physical  phenomenon. 

A  box  is  shown  in  Fig.  2.8(a)^^'^'^  whose  sides  are  of  unit  length.  Next  to  the 
box  is  a  shape  called  a  "generator,"  composed  of  line  segments  scaled  to  1/4  of  the 
length  of  a  side  of  the  box.  Each  side  of  the  square  is  replaced  with  this  "scaled" 
generator.  The  result  is  shown  in  Fig.  2.8(b).  Again,  the  generator  is  scaled  down  to 
the  length  of  each  straight  line  segment  of  this  new  object.  The  straight  line  segments 
of  Fig.  2.8(b)  are  replaced  with  the  scaled  generator.  A  portion  of  the  resulting  shape 
is  Fig  2.8(c).  This  process  is  continued,  generating  an  object  which  is  said  to  be  scale 
invariant  and  self-similar.  A  scaling  fractal  looks  geometrically  the  same  everywhere 
and  on  all  scales. 
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Fractal  on  bounded 
range 


logL 


logE 


Fig.  2.7  Physical  fractals  have  scale  dependence  over  a  bounded  range  of  scales. 
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(a) 


(b) 


(c) 


each  side  of  the  square 
is  replaced  by  the  scaled 
generator. 


again,  the  generator  is 
scaled  to  replace  each 
of  the  segments  above. 


Fig.  2.8  The  stepwise  construction  of  a  scaling  fractal. 
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The  fractal  dimension  of  this  object,  however,  can  be  computed  in  another, 
quite  different  fashion.   Scaling  fractals  obey  the  following: 

Nr°  =  1,  (2.18) 

where  N  is  the  number  of  elements  in  the  generator, 

r  is  the  scale  factor  of  an  element, 
D  is  the  similarity  or  fractal  dimension. 
For  the  plot  of  Fig.  2.8  called  the  quadratic  Koch  curve,  N  =  8,  r=  1/4,  Thus, 
the  fractal  dimension  D  =  log8/log4  =  3/2. 

Fracture  Surface  Analysis  by  Fractal  Geometry 

Mecholsky  et  al.^''''^^  experimentally  determined  a  relationship  between  fracture 
toughness  (K^)  and  fractal  dimension  (D).   They  found  that 

Kc  =  A(D-1)1^2  ^  AD*^^2^  (2.19) 

where  D*  is  the  fractional  part  of  the  fractal  dimension, 

Kc  =  Kic  -  Kq,  where  Kg  is  the  toughness  of  the  material  for  a  smooth 

(Euclidean)  fracture  surface,  and 
A  is  a  constant. 

Thus,  as  fractal  dimension  increases,  fracture  toughness  increases.     In  some 
papers  (D-1)  is  referred  to  as  D*.   Thus  the  above  equation  can  be  rewritten  as 

Ke    =AD*^'^,  (2.20) 

where  A  is  a  family  parameter.  The  family  parameter.  A,  identifies  a  line  in  the  K^ 
vs.  D*  plane.  Within  a  family,  an  increase  in  fractal  dimension  corresponds  to  an 
increase  in  fracture  toughness. 

Since  D*  is  dimensionless,  a  dimensional  analysis  of  this  equation  requires  that 
A  have  the  dimensions  of  toughness.   Mecholsky  et.  al.^^  proposed  that  A  is  a  product 
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of  a  characteristic  length  and  Young's  modulus,  so  that  A  =  E(ao)^^^.   Then  the  above 
equation  can  be  written  as: 

Kc  =  E  (  ao  D*  )l^^  (2.21) 

where  K^  is  the  fracture  toughness, 

E  is  the  Young's  modulus, 
D*  is  the  fractal  dimension,  and 
a.Q  is  the  characteristic  length. 
Recalling  that     Kg  =  VeG  =  ^2yE    for  a  plane  stress  condition,  a^  can  be 

obtained  by  the  following  equation, 

ao  =  ^^.  (2.22) 

0     ED* 

In  order  to  prevent  getting  an  infinite  values  of  ao  calculated  from  the  above 
equation  if  D*  is  approaching  zero.  A  modified  equation  from  a  different  approach  to 
find  ao  is  expressed  as: 

y  =  Yo  +  aoEhkiD*/2,  (2.23) 

where  y  is  the  fracture  energy, 

Yo  is  the  surface  energy  for  an  Euclidean  fracture  surface  (D*  =  0), 

ao  is  a  characteristic  (atomic)  length,  and 

Ejjj  is  the  Young's  Modulus  of  the  hkl  fracture  plane. 

Using  experimentally  obtained  values  of  y,  E  and  D=^,  yo  and  ao  can  be 
calculated.  The  exact  interpretation  is  not  clear  at  this  time;  however,  it  is  reasonable 
to  assume  that  it  is  related  to  the  characteristic  length  of  the  generator.  In  terms  of 
fracture,  this  would  correspond  to  a  characteristic  bond  length  which  is  involved  with 
initial  fracture  at  the  crack  tip. 

Mecholsky  et  al.^  reported  that  the  outer  mirror  constant,  i.e.,  aitx}^'^  where  t2 
is  as  noted  in  Fig.  2.2,  is  related  to  Kc  and  independently  to  E.  Kirchner49  shows  a 
similar  relationship  with  the  crack-branching  constant  and  used  this  information  to 
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support  his  idea  of  a  constant  strain  intensity  at  branching.  Tsai  and  Mecholsky^O  used 
a  fracture  energy  concept  to  prove  that  the  mirror-hackle  boundary  forms  at  a  constant 
energy  value.  Mecholsky^  1  points  out  that  the  rj/c  =  constant.  Freiman  et  al.52 
showed  that  the  relationship  could  be  represented  as: 

Kbj=(bj)l/2E,  (2.24) 

where  b;  is  a  characteristic  dimension  on  the  atomic  scale,  and 

j  =  1,  2,  3  is  related  to  mirror-mist,  mist-hackle  and  crack  branching 
boundary. 

They  further  showed  that  if  bj  =  a^,  where  bj  corresponds  to  the  mirror-mist 
boundary,  then  the  flaw-to-mirror-size  ratio  is  equal  to  the  fractal  dimensional 
increment: 

c/ri  =  D*.  (2.25) 

This  equation  can  provide  a  link  between  fractography  and  fractal  geometry. 
They  found  equation  (2.25)  in  good  agreement  with  experimental  results.  Thus,  it 
appears  that  the  relationship  between  the  crack  size  and  the  formation  of  the  branching 
boundaries  is  related  to  the  dimensionality  of  the  strucUire.  Presumably,  the 
dimensionality  of  the  structure  is  related  to  the  bonding,  i.e.,  the  strength  and  length  of 
bonds  on  the  atomic  scale. 

Molecular  Dynamics  Simulation 

Molecular  dynamics  simulation  is  a  computer  technique  to  model  material 
structure  at  the  atomic  level  in  order  to  model  or  study  the  material  properties  with  an 
atomic  view  point.53,54  During  the  molecular  dynamics  simulation,  a  small  box  with  a 
limited  number  of  atoms  is  generated.  The  small  box  is  considered  a  simplified  system 
which  is  ideally  a  realistic  subset  of  the  real  system.  This  MD  technique  uses 
Newton's  equations  of  motion  at  constant  acceleration  over  very  short  time  intervals. 
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Interaction  forces  generated  from  the  potential  energy  generated  between  the  atoms  are 
assumed.  The  acceleration  of  particles  (or  atoms)  results  from  the  interaction  of  those 
interatomic  forces.  The  interaction  forces  produce  the  fundamental  momentum  for  the 
particles  inside  the  small  and  idealized  system. 

Molecular  dynamics  simulations  can  be  applied  to  study  the  properties  of 
materials.  During  MD  simulations,  the  position  and  velocity  of  each  atom  can  be 
calculated  to  show  the  dynamic  properties  at  that  instant  of  time.  The  power  of  the 
MD  simulation  is  that  the  rule  of  interaction  can  be  changed  at  will  to  model  any 
complex  sets  of  particles,  and  environments  can  be  simulated  at  conditions  which  are 
not  accessible  experimentally.  The  difficulty  with  the  MD  simulation  is  that  the  rule  of 
interaction,  i.e.,  the  selected  potential,  must  be  carefully  chosen  in  order  to  be  realistic 
or  meaningful. 

Three  major  uses  of  MD  simulation  are  the  development  of  theoretical  models 
for  materials,  the  testing  of  theoretical  models,  and  the  predictions  of  material  behavior 
under  experimentally  inaccessible  or  costly  conditions. 

Material  models  can  be  developed  from  molecular  dynamic  simulations  to 
observe  what  occurs  during  the  simulation  and  be  able  to  set  the  rules  of  a  simulation. 
If  we  want  to  model  the  fracture  of  silicon,  we  do  this  by  observing  individual  atomic 
motions  as  samples  are  fractured  using  identical  interaction  potentials.  We  can  also 
fracture  at  several  strain  rates  to  investigate  the  effects  of  allowing  the  material  to  relax 
by  thermal  vibrations  as  compared  to  not  allowing  the  material  to  relax.  After 
performing  such  simulations,  we  can  develop  a  general  idea  of  how  single  crystal 
silicon  fractures  and  those  ideas  will  be  valuable  to  those  studying  the  strength  of  Si. 
Thus,  molecular  dynamics  simulations  can  be  used  to  devise  models. 

Molecular  dynamics  can  be  also  used  to  test  theoretical  models  using  the  ability 
to  observe  all  modeled  atomic  motions  and  the  ability  to  simulate  a  wide  range  of 
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environments.  If  the  system  behaves  as  predicted  by  a  model  under  a  range  of 
conditions,  then  the  model  can  get  greater  credibility. 

Extreme  conditions  also  can  be  simulated  using  molecular  dynamics.  As  one 
might  simulate  any  environment,  one  needs  only  to  have  a  realistic  simulation  to 
predict  materials  behavior  or  study  possible  modes  of  failure  in  extreme  environments. 
Thus,  as  molecular  dynamics  becomes  better  at  simulating  reality,  one  may  expect  this 
application  of  molecular  dynamics  to  find  regular  use. 

When  creating  a  molecular  dynamics  system,  one  would  try  to  model  the  laws 
of  nature  as  close  to  reality  as  possible  in  order  that  the  processes  observed  in  the  MD 
simulation  can  be  as  realistic  as  possible,  and  that  the  resulting  character  of  the 
simulated  system  corresponds  to  the  character  of  real  materials.  Because  the  laws  of  an 
MD  system  are  executed  by  high  speed  computer  of  limited  computational  power,  there 
must  be  limitations  of  the  resulting  MD  system.  These  lunitations  are  seen  as 
limitations  in  the  spatial  and  temporal  resolution  in  nature,  as  well  as  limitations  on  the 
number  of  interactions  possible  within  such  a  system. 

The  spatial  dimensions  of  a  molecular  dynamics  cell  must  be  within  a  certain 

en 

order  of  magnimde  of  the  computer  precision  of  positions  obtainable.  The  spatial 
resolution  of  the  computer  calculation  is  generally  1  in  10'  or  1  in  10  °,  depending  on 
whether  single  or  double  precision  is  used  in  the  computer  calculations.  Over  long 
calculation  periods,  the  use  of  single  precision  has  resulted  in  a  small  error  in  the 
system,  while  yielding  a  lot  of  time-saving  is  obtained.  Thus,  here  the  choice  of  spatial 
resolution  is  not  a  problem  in  performing  the  MD  simulations. 

The  length  of  time  represented  by  the  simulation  must  be  within  a  certain  order 
of  magnitude  of  the  time  step,  the  basic  unit  of  time  in  a  simulation.  Generally  the 
duration  of  an  experiment  is  between  10,000  to  100,000  time  steps  depending  on 
computational  resources  and  the  purpose  of  the  experiment.  The  temporal  resolution  of 
an  MD  simulation  may  be  a  serious  limitation.    The  reason  for  this  is  that  processes 
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relevant  to  the  structural  evolution  of  materials  have  characteristic  times  spanning  about 
19  orders  of  magnitude,  far  more  than  the  four  to  five  orders  of  magnitude  spanned  by 
MD  simulations.  Generally  a  simulation  which  lasts  for  hundreds  or  thousands  of 
atomic  thermal  vibrations  is  needed  in  order  to  understand  the  properties  of  the  system. 

The  mass  of  the  electron  is  about  one  ten  thousandth  of  those  of  atoms.  The 
motion  of  electrons  cannot  be  modeled  with  the  motions  of  atoms  at  the  same  time 
since  the  characteristic  time  is  so  small  compared  with  that  of  atoms.  Thus  electronic 
motions  cannot  be  used  to  provide  true  interatomic  potentials  of  materials  in  such  a 
simulation  and  the  motions  studied  here  are  limited  to  those  of  atoms  only.  The  chosen 
potential  is  necessary  to  accurately  represent  the  thermal  vibration  motions  of  atoms  in 
order  to  model  the  properties  of  the  simulated  system.  Thus  the  smallest  time  step  used 
must  be  a  small  fraction  of  the  vibration  period.  A  typical  time  step  chosen  in  this 
study  is  0.5  x  10'^^  sec.  One  generally  can  only  run  simulations  for  about  10^  time 
steps  which  results  in  a  limitation  in  the  duration  of  modeled  experiments  to  hundreds 
or  thousands  of  atomic  vibrations. 

Another  limitation  is  that  the  number  of  interactions  between  atoms  cannot  be 
too  big.  The  more  atoms  put  in  the  system,  the  more  computer  CPU  time  will  be 
consumed.  For  example,  if  a  system  with  n  particles  is  created  and  only  one-to-one 
interaction  is  considered,  then  the  total  interactions  between  those  atoms  is  n^.  The 
time  to  calculate  a  system  with  lOn  atoms  is  10^  longer  than  that  needed  to  calculate 
one  with  n  atoms.  Thus  the  simulated  system  must  be  kept  small  and  the  chosen 
interacted  potential  must  be  simple  and  easily  applied. 

Thus  the  differences  between  an  MD  system  and  the  real  universe  are  that  the 
length  of  studied  time  and  the  size  of  the  studied  system.  The  MD  system  has  to  be 
small  because  of  the  limits  on  the  spatial  resolution  of  time  and  number  of  atoms  used 
in  the  system. 


CHAPTER  3 
MOLECULAR  DYNAMICS  TECHNIQUES 


Overview 

A  technique  for  studying  materials  whose  popularity  and  usefulness  has  arisen 
within  the  past  few  decades  is  molecular  dynamics  (MD)  computer  sunulation.  This 
technique  can  solve  equations  of  motion  for  a  system  of  particles  which  may  be  ions 
(ionic  bonding)  or  atoms  (covalent  bonding  or  metallic  bonding).  Each  particle  moves 
according  to  the  forces  on  it  caused  by  all  of  the  other  particles  in  the  system.  The 
inter-particle  forces  are  derived  from  an  assumed  potential  function  which  describes  the 
interactions  between  the  atoms.  Several  types  of  systems  have  been  simulated  by  MD 
with  various  interatomic  potentials.  Those  systems  include  hard  sphere  systems,  soft 
sphere  systems,  ionic  systems  and  Lennard- Jones  systems.  These  systems  are 
generally  simulated  using  simple  pair  interaction  potentials  to  obtain  a  close 
approximation. 

The  covalent  bond  formed  by  two  Si  atoms  sharing  electrons  is  very  localized 
and  directional.  In  ionic  systems.  Coulomb  potentials  which  mclude  short  range 
repulsive  forces  as  well  as  attractive  forces,  are  required.  In  covalent  systems,  bond- 
directionality  should  be  included.  Two  body  interatomic  potentials  are  not  enough  to 
produce  the  diamond  cubic  strucmre.  Since  Si  exhibits  that  structure  as  a  solid,  three 
body  potentials  must  be  introduced  in  order  to  conduct  a  reasonable  simulation.  The 
Stillinger-Weber  potential^"^  has  been  used  for  Si  solid  simulations  using  MD  and  the 
simulated  pair  correlation  function  agrees  very  well  with  the  experimental  data.    It  was 
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34 
also   generally  found   able  to   accurately   simulate   the   elastic   and   thermodynamic 
properties  of  Si.^^'^^   This  potential  has  been  adopted  here. 

Potential  Determination 

Single  crystal  silicon  consists  of  atoms  held  in  place  (a  diamond  structure  which 
is  relatively  anisotropic  from  the  atomic  viewpoint)  by  strong  and  directional  bonds.  It 
seems  reasonable  at  first  sight  that  the  corresponding  potential  O  could  be 
approximated  by  a  combination  of  pair  and  triplet  potentials,  O2  and  O3.  And  these 
two  potentials  can  be  expressed  as  a  function  of  energy,  s,  and  length  units,  5.  Thus, 

cD2(rij)  =  s  f2(rij/5),  (3.1) 

and 

03(r„  rj,  rk)  =  s  U^/b,  rj/5,  rk/5).  (3.2) 

We  use  a  reduced  pair  potential  for  Si  selected  from  the  following  five- 
parameter  family,  a  simplified  Lennard- Jones  potential,  from  the  work  of  Stillinger  and 
Weberi^"^ 

f2(r)   =  A  (B  r'P  -  r-^)exp[(r  -  a)"^],    if  r<a;  (3.3) 

=  0,   if  r>a. 
Thus, 

02(rij)  =  s  f2(rij/§)  (3.4) 

=  8A(B(^)-P  -  (^)-'')exp[(  '-  -  J)-l],    ifr<a; 
80  00 

=   0,   ifr>a. 

This  form  automatically  cuts  off  at  r  =  a  without  discontinuities  in  any  r  derivatives, 

which  is  a  distinct  advantage  in  any  MD  simulation  application.    Stillinger  and  Weber 

use  the  same  cutoff  technique  for  the  three-body  interactions,  f3(ri,rj,rk)  =  hjjk  +  h^k 

-I-  hkij.   The  function  h  is  given  : 
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^2 


h(rij,  %,  Bjik)  =  X  exp[ti  /(rij  -  a)  +  r|  /(rjk  -  a)]  x  (cosBjik  +  1/3)  ,   if  r<a; 

(3.5) 

=  0,  ifr>a. 
where  Gjik  is  the  angle  between  Tj  and  r^.    The  function  h  exists  when  both  ly  and  r^ 
are  less  than  the  previously  introduced  cutoff  a. 

The  parameters  which  best  satisfy  the  radial  distribution  function  were 
determined  by  Stillinger  and  Weber  to  be:  A  =  7.049556277,  B  =  0.6022245584,  p  = 
4,  q  =  0,  a  =  1.8,  1  =  21.0,  h  =  1.20,  d  =  0.20951  nm,  and  s  =  50  kcal/mole  = 
3.4723x10"^^  erg/atom  pair. 

Although  the  Stillinger-Weber  potential  has  been  successfully  used  to  model  the 
strucmre  of  silicon,  it  does  not  provide  the  expected  repulsive  forces  following  covalent 
bond  breakage  which  result  from  ion  formation  during  the  fracture  process.  The  Si-Si 
covalent  bond  is  formed  by  an  attractive  force  even  if  the  distance  between  the  two 
atoms  is  larger  than  the  equilibrium  position  as  long  as  it  is  within  the  cutoff  distance. 
However,  during  fracture  this  condition  is  not  realistic  and  makes  the  Si-Si  bond 
difficult  to  break.  An  unrealistically  long  fiber-like  strucmre  is  observed  to  result  if  the 
unmodified  Stillinger-Weber  potential  is  used.  An  improvement  developed  here 
assumes  that  following  bond  breakage  each  Si  atom  will  gain  one  positive  charge  for 
each  lost  neighbor.  Si  will  have  two  positive  charges  if  it  loses  two  of  its  four 
neighbors  and  four  positive  charges  if  all  four  neighbors  are  displaced  beyond  the 
cutoff  distance. 

A  modified  Born-Mayer  potential^^'^^  is  used  for  modeling  the  resulting  ionic 
forces.  Consequently  this  approach  assumes  ionic  repulsive  forces  when  two  Si  atoms 
both  with  positive  charges  fall  inside  the  Born-Mayer  cutoff  distance.  This  potential 
has  the  following  form: 
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r^ij  =  Aij-exp(— )  +  (^-erfc(-f ),  (3.6) 

where  erfc()  is  the  complement  error  function, 

s  is  a  constant, 

L  is  a  linear  dimension  of  the  molecular  dynamics  cell, 

q  is  the  charge, 

r  is  the  distance  between  the  i  and  j  atoms,  and 

p  is  a  hardness  parameter. 
After  combining  the  Stillinger-Weber  potential  and  the  modified  Born-Mayer 
potential,  the  fiber-like  structure  disappears  and  potentially  more  realistic  fracture 
surfaces  can  be  obtained. 

Fundamentals  About  the  MP  Simulation 

How  the  MD  simulations  are  performed  will  be  explained  now.  First,  we 
discuss  how  the  system  was  set  up,  which  includes  the  initial  particle  positions, 
velocities  and  the  potential  described  in  the  previous  section.  Then  we  discuss  how  the 
system  is  operated,  which  includes  the  application  of  periodic  boundary  conditions,  the 
force  summation  on  the  particles,  control  of  the  velocities  of  particles,  and  control  of 
thermodynamic  values.  Finally  the  discussion  of  how  the  data  from  the  simulation 
system  is  collected  will  be  given.  A  flow  chart  in  Fig.  3.1  is  given  to  illustrate  how 
the  MD  sunulation  is  operated. 

Initial  Conditions 

Initial  particle  positions  are  located  at  specific  positions  since  a  single  crystal 
structure  is  desired. 


37 


Initial  Velocity: 


temperature 


kinetic 
energy 


random 

distribution 

function 


velocity  (V) 


Interactions  between  atoms: 


apply  strain 


position 
of  atoms 


+  (V  +  A  V)  X  A  t 


new  position 
of  atoms 


1 


acceleration 


AV 


force 

between 

atoms 


ax  At 


Fig.  3.1  A  flow  chart  for  the  MD  simulation. 
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Initial  velocities  are  given  using  the  Poisson  velocity  distribution.  The  velocity 
distribution  is  generally  selected  to  correspond  to  a  temperature  at  which  the  crystal  is 
stable.  Here  it  corresponds  to  room  temperature  or  to  the  temperature  of  the 
simulation. 

Random  velocities  are  assigned  to  particles  using  the  Poisson  velocity 
distribution.      The   formula  giving   the  Maxwellian  velocity   distribution  in  three 

CO 

dimensions  are 


Vxi  = 

^^^^  X  In(Randl) 

Vy.= 

^'^^^  X  In(Randl) 

Vzi  = 

^^^^xln(Rand2) 

X  cos(27r  X  Rand3) 


X  sin(27i  X  RandB) 


X  sin(27xx  Rand4) 


(3.7) 


(3.8) 


(3.9) 


where 


Vxi,  Vyj,  and  V^i  are  the  components  of  the  velocities  of  particle  i, 

Randl,  Rand2,  Rand3,  and  Rand4  are  random  numbers  evenly 

distributed  between  0  and  1 , 

mj  is  the  mass  of  particles  i, 

T  is  the  temperature,  and 

kg  is  the  Boltzmann's  constant. 
Particle  velocities  are  applied  by  putting  the  velocities  obtained  above  into 
positions  before  and  after  the  current  time  step,  which  is  the  basic  unit  of  time  during  a 
simulation.  Thus  the  random  velocities  are  effectively  converted  into  average  velocities 
during  the  previous  time  step,  and  particles  positions  and  velocities  are  given  by  the 
difference  of  two  sets  of  positions  for  each  particles.  The  reason  for  this  will  become 
clear  in  the  discussion  of  Verlet's  algorithm. 
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Tnter actions  Between  Atoms 


Interaction  forces  are  tabulated  for  all  types  of  interactions  possible  over  a  range 
of  interparticle  spacing.  The  interatomic  potential  assumed  in  this  study  is  the 
Stillinger-Weber  potential  and  the  modified  Born-Mayer  potential  as  illustrated  above. 
The  derivative  of  these  potentials  is  the  force  as  a  function  of  r.  The  two-body 
potential  is 

02(rij)  =  s  f2(rij/§) 

=  8A  (B  (|)-P  -  (^)-'l)exp[(  ^  -  ^y'],  if  r<a; 
§  6  0       0 

=  0,  if  r>a. 
Thus,  the  two-body  force  is 
5O2 


f2(rii) 


dr 


eA 
=  —  exp 
S 


r-a 


y' 


(3.10) 


(Bp(^)-P-l-(^)-^-^)  + 
o  0 


B 


v5. 


^-q^ 


52 


(r-a)^ 


,  if  r<a; 
=  0,  if  r>a. 
It  is  computer-time-effective  to  tabulate  values  of  this  function  as  a  function  of 
r^  and  divide  by  r  so  that  one  may  use  the  square  of  the  interaction  distance  in  finding 
interparticle  forces  in  the  table.  One  may  find  the  Cartesian  components  of  the  force 
between  two  particles  by  simply  multiplying  the  appropriate  tabulated  value  by  the 
difference  in  Cartesian  coordinates  of  the  two  particles. 
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The  energy  of  potential  energies  for  interactions  between  particles  are  also 
tabulated  in  terms  of  r^.    This  is  used  for  the  calculation  of  potential  energies  for  the 
system. 

The  three-body  force  is  considered  in  the  same  way. 
(D3(ri,  rj,  ry.)  =  s  fsCq/S  ,  rj/6  ,  rk/6). 
f3(ri,rj,rk)  =  hji,,  +  hy^  +  hkij- 
The  function  h  is  given  : 

hjik  =  h(rij,  rjk,  Ojik)  (3.11) 

=  X  exp[ri  /(rij  -  a)  +  r)  /(rik  -  a)]  x  (cosBjik  +  1/3)^,   if  r<a; 
=  0,   if r>a. 
The  following  derivation  related  to  the  three-body  potential  will  only  use  hji^  for 
convenience.   Thus  the  three  body  force  related  to  hji^  is  as  follows: 

j^    -  SO    ^   ISO  .^  i^v 

VO  =  gradO  =  e^  --  +  eg  -— -  (3 .  12) 

or  r  99 

,^  ao     ^1  50,  .^  i-,s 

F  =  -VO  =  -(er  — +  ee-— )  (3.13) 

dr  r  50 

(A)  At  first,  we  differentiate  r^  only. 

F  =  Fr+Fe=-VO  (3.14) 

At  Cj.  direction, 

K  =-^  =  -^cxp[jin-a)-'+y(r2-a)-'](cosQ  +  \f  (3.15) 

^         dri      {ri-ay  ^ 

At  §9  direction, 

FQ=-l^  =  ^exp[Y(ri-a)-l+Y(r2-a)-^]sine(cose  +  ^)  (3.16) 

rj  99       rj  3 

Now  transform  the  polarized  coordinates  into  Cartesian  coordinates. 

Fr+Fe=Fx+Fy+F,  (3.17) 
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Consider  F^  first.    At  first  we  have  to  find  the  transformation  of  u^  mto  Cartesian 
coordinate. 

ui=U]xT  +  uiyj  +  uizk  (3-18) 

Friii=Fr(uiJ  +  uiyj  +  u,2lc)  (3.19) 

Now  consider  Fq.   According  to  Fig.  3-1, 

60  =ui  XU3  =U]  x(ui  XU2)  (3.20) 

and  because  A  x  (  B  x  C  )  =  (A-C)xB-(A-B)xC,  thus 

ee  =  uix(iii  XU2)  (3-21) 

=  (iil-U2)xui-(ui-ui)u2 

=  cosBuj  -U2 
Fgee  =  Fq (cos Bill  -112)  (3.22) 

=  Fe(coseuix  -  U2x)i  +  Fe(coseuiy  -  U2y) j  +  F0(coseuiz  -  U2z)k 

Thus, 

Fr+Fe  =  Fx+Fy+F,  (3.23) 

=      [FrUix  +  Fe(coseuix-U2x)]i 
f[FrUiy  +  Fe(coseuiy-U2y)]j 
f[FrUiz  +  F0(coseu,z  -  U2z)]k 
(B)  We  can  consider  the  r2  term  in  the  same  way  as  (A). 

MP  Simulation  Procedures 

The  explanation  of  how  the  simulations  are  performed  will  be  given  here.53 
The  concept  of  time  step,  the  interaction  between  particles,  the  method  used  to  update 
particle  positions,  the  application  of  periodical  boundary  conditions,  and  the  control  of 
system  temperature  will  be  discussed  as  follows. 
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The  Concept  of  Time  Step 

The  time  step  is  the  basic  unit  of  time  during  a  molecular  dynamics  simulation. 
As  there  are  many  atoms  interacting  in  the  simulated  system,  the  motions  of  particles 
cannot  be  solved  analytically.  Thus  one  can  calculate  the  motions  of  atoms  using 
Newton's  laws  of  motion  where  a  time  step  dt  is  used  in  which  the  acceleration  is  kept 
constant.  A  careful  choice  of  time  step  is  very  important.  If  the  time  step  is  chosen  to 
be  too  small,  a  waste  of  computer  running  time  results.  If  the  time  step  is  chosen  to  be 
to  big,  nonlinear  effects  have  to  be  considered  and  the  atoms  will  move  too  far  away  to 
reach  unrealistic  positions  before  the  next  time  step  is  calculated.  Following  the  work 
of  Soules,  Ochoa,  and  Swiler,  the  time  step  used  in  our  simulation  is  0.5  x  10'^^ 
second,  and  is  about  1/40  of  an  atomic  vibration  period.  A  test  to  decide  if  a  time  step 
is  of  sufficient  length  is  to  perform  a  simulation  over  about  a  thousand  steps  and  check 
to  see  if  the  system  energy  remains  fairly  constant.  If  not,  the  time  step  is  too  long.  A 
time  step  of  0.5  x  10'^^  seconds  was  found  to  be  of  acceptable  length  and  was  applied 
here. 

The  Method  to  Update  Particle  Positions 

It  is  very  important  to  handle  particle  interactions  in  a  time-saving  way. 
Particle  interactions  are  allowed  only  when  particles  are  within  a  certain  distance  for 
others.  Beyond  the  distance,  the  effect  of  interaction  is  either  none  or  negligible. 
Cutoff  distance  is  the  distance  in  which  interactions  are  allowed.  Thus  a  neighbors-list 
for  all  particles  which  lists  all  the  neighbors  which  may  be  considered  for  the  next  5  or 
10  time  steps  will  be  made  and  updated  at  every  5  or  10  time  steps  depending  on  the 
pulling  rate.  The  fixed  distance  for  the  neighbors-list  is  chosen  to  be  the  sum  of  the 
cutoff  distance  and  the  distance  a  particle  may  move  in  the  next  5  or  10  time  steps. 
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Only  the  neighbors  Usted  on  the  neighbors-list  will  be  considered  when  calculating 

interactions. 

The  squares  of  the  distances  between  particles  are  calculated  using  the 
Pythagorean  theorem  for  time-saving  reasons.  Square  roots  of  these  values  are  not 
taken  when  calculating  interaction  forces  because  square  root  operations  are  time- 
consuming  and  the  force  table  is  tabulated  in  terms  of  distance  squared.  Square  roots 
of  the  distances  squared  are  only  required  when  the  neighbors-list  are  updated  and 
when  pair  correlation  functions  and  bond  angle  distribution  functions  are  accumulated. 

Periodic  Boundary  Condition 

Periodic  boundary  conditions  are  applied  when  finding  inter-particle  distances. 
If  periodic  boundary  conditions  are  applied,  particles  which  leave  the  primary  cell  after 
a  time  step  will  be  translated  back  into  it  by  subtracting  the  appropriate  cell  dimension. 
Thus  particles  which  leave  through  one  side  of  the  system  cell  will  appear  through  the 
opposite  side  of  the  cell. 

Thus    if   the    difference    in   the    coordinate    is    greater    than    one-half   the 

corresponding  dimension  of  the  system  cell,  the  corresponding  image  of  the  particle 

located  in  the  adjacent  cell  will  be  closer  to  the  central  particle  than  the  particle  in  the 

primary  cell.    The  difference  in  the  coordinate  is  then  adjusted  to  reflect  this  by  either 

adding  or  subtracting  the  corresponding  dimension  of  the  cell  to  minimize  the  absolute 

values  of  the  difference.   It  is  the  resulting  difference  which  is  used  in  the  Pythagorean 

relation  to  find  the  square  of  distance  between  the  particles.    The  algorithm  by  which 

this  is  performed  is  as  follows: 

dijk=Xij-Xik  (3.24) 

d'^,.d,,-s,.™nc(^)  0.23) 
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4=V'+'^2jk'+d'3jk'  (3-26) 

where  i  =  1,  2,  3  denote  the  Cartesian  dimension  of  space, 

Xjj  is  the  i  Cartesian  coordinate  of  atom  j, 
Xy^  is  the  i  Cartesian  coordinate  of  atom  k, 
djjk  is  the  i  coordinate  difference  between  atoms  j  and  k, 

d'..,   is  the  i  coordinate  difference  between  atoms  j  and  k  after  correction 

ijk 

for  periodic  boundary  conditions, 

Sj  is  the  i  dimension  of  the  system  cell, 

trunc( )  is  a  function  which  returns  the  integer  closest  to  the  argument 

between  the  argument  and  zero,  and 

rjj  is  the  distance  between  atoms  i  and  j . 


Interactions  Due  to  Two-Body  Potential 

Two-body  interactions  between  particles  are  handled  in  a  four  step  process. 53 
The  square  of  the  interparticle  distance  is  found  for  particles  located  in  the  neighbors 
lists.  The  interparticle  force  is  then  found  by  looking  up  the  appropriate  value  in  the 
force  table.  The  components  of  the  force  are  found  by  multiplying  the  value  tabulated 
in  the  force  table  by  the  difference  in  the  coordinates  for  each  dimension  (dyk).  The 
components  of  force  are  then  summed  for  each  particle  for  all  interactions. 

Interactions  Due  to  Three-Rndv  Potential 

Three-body  interactions  between  particles  are  handled  in  a  four  step  process. 
The  square  of  the  interparticle  distance  is  found  for  particles  located  in  the  neighbors 
lists.  Two  atoms,  j  and  k,  within  the  cutoff  potential  range  are  found  from  the 
neighbor  list  for  the  center  atom,  i.    The  angle  Gjjk  between  the  center  atom  i  and  the 
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other  two  atoms  j  and  k  were  calculated.  The  interparticle  force  is  then  found  by  the 
derivative  of  the  three-body  potential.  The  radial  force  and  the  tangential  force  are 
calculated  separately.  The  components  of  the  force  in  each  dimension  i,  j,  and  k  are 
found  by  multiplying  the  previous  calculated  tangential  and  radial  forces  with  respect  to 
the  related  vectors  in  the  coordinates.  The  components  of  force  are  then  summed  for 
each  particle  for  all  interactions. 

Updating  of  the  Atom  Positions 

After  the  forces  on  each  atom  are  accumulated  from  either  the  two-body  or  the 
three-body  potential,  the  change  of  positions  of  particles  are  updated  using  Verlet's  or 
Gear's  algorithm  at  each  time  step.  The  details  of  the  Verlet's  algorithm^^  and  Gear's 
algorithm^^  will  be  illustrated  as  follows. 

Verlet's  algorithm 

The  advantages  of  Verlet's  algorithm59  are  speed  and  small  error.  The 
derivation  for  this  algorithm  is  as  follows:  According  to  Taylor  series  expansion,  the 
position  at  time  (t  +  At)  calculated  from  time  t  is 


Xj(t  +  At)  =  X(t)  + 


dXi(t) 


At.i 

dt    ;      2 


r^2 


d"Xi(t) 
dt2 


At2+1 
3! 


dt^ 


At^+OCAt"^); 


while  the  position  at  (t  -  At)  calculated  from  time  t  is 


(3.27) 


Xi(t-At)  =  Xi(t)-l^]At  +  l 


(a! 


d^Xi(t) 
dt2 


At^-l 
3! 


^dbq(t)^ 
dt- 


At^+OCAt"^). 


(3.28) 
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Adding  these  two  expansions,  we  can  get  the  new  position: 

At2+0(At4),  (3.29) 


Xi(t  +At)  =2Xi(t)-Xi(t-At)  + 


^d^XiCt)^ 


V     dt^     J 


^hej-e  i  =  1 ,  2,  3  denote  the  Cartesian  dimensions  of  space, 

t  is  a  particular  time, 

Xi(t)  is  the  i  Cartesian  coordinate  at  time  t, 
Xi  (t  -  At)  is  the  i  Cartesian  coordinate  at  time  t  -  At ; 
At  is  the  length  of  a  time  step,  and 
0(At'^)is  the  remainder  terms. 

Gear's  algorithm 

Gear's  algorithmic  is  another  choice  to  calculate  the  change  of  positions  of 
particles  due  to  two-body  interaction.  The  advantage  of  Gear's  algorithm  is  that  it 
gives  less  energy  fluctuation;  while  the  big  disadvantage  is  that  it  is  very  time- 
consuming.  Gear's  algorithm  predicts  molecular  positions  Xj  at  time  (t  -I-  At)  using  a 
fifth-order  Taylor  series  based  on  positions  and  their  derivatives  at  time  t.  Thus,  the 
derivatives  Xj^^^,  Xi^"\  Xj^"^),  Xj^^^l  x/^>  are  needed  at  each  step;  these  are  also 
predicted  at  time  (t   -h  At)  by  applying  Taylor  expansions  at  time  At. 

Xi(.  +  At)  =  Xi(t)  +  Xiffl(t)At  +  Xi(«)(t)^  +  Xi(-)(t)*f^ 

+x/")(t)^  +  XiW(t)^  (3.30) 

Xi«  (t  +  At)  =  Xi«  (t)  +  Xi(")  (t)  At  +  Xi(-)  (t)  ^  +  X/-)  (t)  ^ 

+x/^)(t)^  (3-31) 

4! 
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9  a 

Xi(")  (t  +  At)  =  Xi(")  (t)  +  Xi("^)  (t)  At  +  X/*^)  (t)  ^  +  Xi(^)  (t)  ^    (3 .  32) 
x/i")  (t  +  At)  =  Xi^"^)  (t)  +  Xi('^)  (t)  At  +  Xi(^>  (t)  ^  (3.33) 


t 
2!     ■  "'      '"     3! 

2! 

Xi^^^^  (t  +  At)  =  Xi^'^)  (t)  +  Xi^^)  (t)  At  (3 .34) 

x/^)(t  +  At)  =  Xi(^)(t)  (3.35) 

In  order  to  correct  the  predicted  positions  and  their  derivatives  using  the 

discrepancy  between  the  predicted  acceleration  and  that  given  by  the  evaluated  force. 

The  force  at  (t  +  At)  obtained  from  Newton's  second  law  can  be  used  to  determine  the 

acceleration  Xi^"^(t  +  At).     The  difference  between  the  predicted  accelerations  and 

evaluated  accelerations  is  then  formed, 

AXi^")  =Xi(")(t  +  At)-Xi(")P(t  +  At).  (3.36) 

In  gear's  algorithm  for  second-order  differential  equations,  this  difference  term 

is  used  to  correct  all  predicted  positions  and  their  derivatives.   Thus, 

Xi=xf+aoAR2  (3-37) 

Xi«At  =  Xi('^PAt  +  aiAR2  (3-38) 

(ii)(At)^^      (ii)p(At)^_^        j^  (3.39) 

'         2\  2! 

(iii)  (Mf_  ^  x.(iiOP  (^  +  a3AR2  (3-40) 

1  31  1  3! 

(,^)  (Mf_  ^     (iv)P  (At)^  ^  ^    .  R  (3.41) 

1  4!  '  4!  ^      -" 

(,)(At)^^      (v)p(At)^^^^R  (3.42) 

5!  '  5!  ' 

2 

Where  AR2=AXi^")^. 

2! 

The  parameter  aj  promotes  numerical  stability  of  the  algorithm.      The  a^ 
depends  on  the  order  of  the  differential  equations  to  be  solved  and  on  the  order  of  the 
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Taylor  series  predictor.  Gear  determined  their  values  by  applying  each  algorithm  to 
linear  differential  equations  and  analyzed  the  resulting  stability  matrices.  For  a  q  order 
predictor,  the  values  of  the  a^  were  chosen  to  make  the  local  truncation  error  of 
0(At'l''"^).  Values  of  the  aj  for  fifth-order  predictors  are  Uq  =  3/16,  a^  =  251/360, 
a2  =  1,  a3  =  11/18,  a4  =  1/6,  a5  =  1/60. 

Temperature  Calculation 

The  temperature  of  the  system  is  controlled  by  assuming  that  the  system 
behaves  as  a  classical  statistical  system,  with  average  particle  kinetic  energies  given  by 
the  equipartition  theorem.  The  temperature  of  the  system  is  then  calculated  by  solving 
the  expression  E^  =  3kT/2  for  T.  The  temperature  is  then  adjusted  by  scaling  the 
velocities  of  all  particles  to  give  the  proper  average  kinetic  energy. 

M.D.  Video  Presentation 

Information  is  collected  from  the  simulation.  Any  aspect  of  the  simulation  can 
be  recorded,  but  there  are  a  few  which  are  commonly  recorded.  These  are  periodic 
snapshots  of  the  structures,  instantaneous  or  average  thermodynamic  properties,  pair 
correlation  functions,  and  bond  angle  distribution  functions. 

Snapshots  of  the  strucmre  of  the  system  are  made  by  periodically  saving  all 
positions  of  particles  of  the  system  in  a  data  file.  These  snapshots  of  the  structures 
may  then  later  be  used  to  perform  a  free  volume  sphere  analysis  on  the  structure  to 
make  graphic  representations  of  the  strucmre. 

Thermodynamic  properties,  such  as  system  temperamre,  energy,  and  pressure 
are  calculated.  Average  of  these  values  may  also  be  kept  in  a  record  file.  Thus  either 
instantaneous  thermodynamic  properties  or  average  properties  may  be  made. 
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Pntential  Enargv  nf  the  System 

The  potential  energy  for  a  system  is  calculated  in  a  way  similar  to  the 
calculation  of  forces  on  individual  particles.  When  the  neighbors-list  is  updated, 
interaction  potentials  are  found  using  the  potential  energy  table.  The  potential  energies 
for  all  interactions  are  then  summed  to  find  the  potential  energy  for  the  system.  The 
system  potential  energy  is  added  to  the  system  kinetic  energy  as  previously  obtained  to 
find  the  system  internal  energy. 

System  Pressure  Calculation 

The  system  pressure  or  stress  tensor  is  calculated  using  the  Virial  of  Clausius, 
which  has  the  form 


P  =  p-k-T 


^      k-n-T 

1  :  ^y^j 


•  Sfe  VO(r)) 


When  converted  for  use  in  simulations,  it  takes  the  form 


10^       ^  ^ 


3V 


mjv^+S^j  Fij 


^  J  J 


where  P  is  the  system  pressure  in  Kbar, 

V  is  the  volume  of  the  cell  in  A  , 

V  is  the  velocity  of  the  particle  in  A/lO'^^s, 
F  is  the  interatomic  force  in  10"     erg/A, 

r  is  the  interatomic  separation  in  A,  and 
mi  is  the  mass  of  the  particle  in  10'     g. 


(3.43) 


(3.44) 
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Pair  Correlation  Function 

Pair  correlation  fiinction  distributions  are  accumulated  when  updating  neighbors 
lists.  These  are  done  by  incrementing  arrays  existing  for  each  type  of  interaction  in 
which  the  elements  of  the  arrays  correspond  to  a  segment  of  distance  every  time  a  pair 
with  some  separation  is  found.  The  contents  of  the  arrays  are  then  converted  into 
distributions  by  normalizing  using  the  function 

G(r).   ^-""''^-^   .  (3.45) 

Nij-47rr^-Ar 

where  G(r)  is  the  pair  correlation  function, 

r  is  the  separation  distance, 
Ar  is  the  separation  distance, 
V  is  the  volume  of  the  cell, 

Ny  is  the  number  of  pairs  found  between  the  distances  r  and 
Ar  from  the  central  atom, 
Nij=m-ni-nj,  (3-46) 

where  m  is  the  number  of  accumulations, 

Hi  is  the  number  of  atoms  of  type  i, 
nj  is  the  number  of  atoms  of  type  j,  and 

47Tr-Ar-^  is  recognized  as  the  volume  of  the  shell  in  which  the  non- 
central  atom  may  be  located  about  the  central  atom. 
This  value  is  output  to  a  normalized  distribution  to  a  data  file. 


51 


Bond  Angle  Distribution 

Bond  angle  distributions  are  found  by  saving  the  atom  identification  numbers  of 
the  coordinating  atoms  about  central  atoms  during  updates  of  the  neighbors  lists,  then 
applying  the  definition  of  a  dot  product  to  find  bond  angles.  The  formula  for  doing 
this  is 


G  =  cos 


XjjXik+yijyik+ZijZjk 
"■ij-fik 


(3.47) 


where  6  is  the  bond  angle, 

i  denotes  the  apex  atom, 

j,  k  denotes  the  end  atoms, 

Xjj,  yij,  Zy  are  the  vector  components  between  i  and  j, 

Xjk,  y-±,  z±  are  the  vector  components  between  i  and  k, 

r;;  is  the  distance  between  atoms  i  and  j,  and 

rjk  is  the  distance  between  atoms  i  and  k. 

Periodic  boundary  conditions  were  used.  An  energy  distribution  function  is 
used  to  assign  the  energy  (velocity)  of  each  atom  at  the  beginning.  The  temperature  of 
the  sample  can  be  varied  by  scaling  the  velocity  of  each  particle  with  an  appropriate 
factor  greater  than  one  and  letting  the  total  kinetic  energy  of  the  particles  be  equal  to 
the  enthalpy  of  the  desired  temperamre.  The  temperamre  of  the  sample  was  kept 
constant  by  scaling  the  velocity  every  time  step  to  avoid  heating  (or  cooling)  effects. 


52 


Tntroduction  of  Strain 

The  determination  of  the  stress-strain  behavior  will  provide  an  insight  into 
material  behavior  for  different  loading  conditions.  Uniaxial  strains  will  be  applied  by 
scaling  the  x-component  of  all  particle  positions  using  the  relation 

Xi  =  Xoi(l+Li/SJ,  (3-48) 

where  Xj  and  x^i  are  the  scaled  and  initial  x-components  of  the  positions  of 

particle  i,  respectively; 

Sx  is  the  x-dimension  of  the  cell  at  the  time  of  scaling,  and 
Li  is  the  desired  incremental  expansion  factor  of  the  cell  m  the  x- 
direction. 
Various  strain  rates  can  be  applied  by  varying  both  the  value  of  Lj  and  the 
frequency  with  which  the  positions  are  scaled.  ^°     In  the  calculation,  the  y  and  z 
directions  will  be  kept  constant,  the  change  in  volume  is  directly  proportional  to  the 
change  in  the  x  direction,  Sx. 


CHAPTER  4 
EXPERIMENTAL  AND  SIMULATION  PROCEDURES 


Experimental  Procedure 

In  this  section,  the  experimental  procedure  is  presented.  Sample  preparation, 
fracture  surface  analysis,  toughness  measurements,  and  fractal  dimension  determination 
will  be  given  in  detail. 

Sample  Preparation 

Single  crystal  silicon  was  provided  by  AT&T  and  IBM.  The  Laue  back 
scattering  method  was  used  to  determine  the  orientation  of  the  crystal.  Several  low 
index  planes  were  chosen  to  be  the  desired  fracture  planes.  They  are  the  {100},  {110} 
and  {111}  planes.  After  the  orientation  was  determined,  single  crystal  silicon  was  cut 
into  flexural  bars  with  the  desired  orientations.     The  surface  of  the  specimen  was 

polished  to  a  Ip-m  finish. 

A  controlled  flaw  was  introduced  in  the  desired  crystalline  plane  using  a 
Vickers  diamond  pyramid  indentation  at  the  center  of  the  tensile  surface  on  flexure 
bars.  Thus,  the  flaw  was  oriented  perpendicular  to  the  longitudinal  axis  of  the 
specimen.  For  the  {110}  orientation,  0.7,  0.9,  1.3,  1.5,  3,  4,  and  5  Kgw  indentation 
loads  were  applied  to  produce  a  controlled  flaw  in  each  bar.  For  the  {100}  and  {111} 
orientations,  1  and  2  Kgw  indentation  loads  were  applied.  Figure  4.1  shows  a 
schematic  arrangement  to  demonstrate  the  location  of  the  indentation  on  the  tensile 
surface.  Flexural  testing  was  performed  on  an  Instron  testing  machine  in  air  at  room 
temperature.    Loads  to  failure  were  recorded.    Fracture  surfaces  were  examined  to 
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Figure  4. 1       Three  orientations  are  chosen  for  comparison. 
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insure  the  samples  had  failed  at  the  indent  site.    The  stress  to  failure  was  calculated 
from  the  dimension  of  each  sample: 

^  -i^k  (4.1) 

^  '  2bh2  ' 
where  a^  is  the  frachire  strength, 

P  is  the  load  to  failure, 
L  is  the  load  span, 
b  is  the  width  of  the  bar,  and 
h  is  the  height  (thickness)  of  the  bar. 

Fracture  Surface  Analysis 

The  procedures  for  fracture  surface  analysis  have  been  well  described  in  the 
literature.6  The  size  of  the  critical  flaw  was  determined  by  modeling  the  flaw  as  an 
idealized  elliptical  crack  of  depth  a  and  half  width  b.  The  equivalent  semi-circular 
crack  radius,  c,  is  determined  by  c  =  Vab  and  is  shown  in  Fig.  4.2.30 

The  fracture  surfaces  were  observed  using  a  light  microscope  to  locate  the 
flaws.  Failure  origins  were  located  by  observing  particular  fracture  markings  on  the 
fracture  surfaces. 61  Crack  branching  and  mirror  features  were  observed  and  were 
used  to  locate  the  flaw. 

Toughness  Measurements 

From  the  measurement  of  the  critical  flaw  size,  c,  and  fracture  strength,  <5^,  the 
critical  stress  intensity  factor,  K^,  for  those  with  residual  stress  caused  from  indentation 
can  be  determined  using 

K,  =  1.65a,c^/2  (4.2) 
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Fig  4  2  Micro-indentation  cracked  bar  under  bending,    (a)  A  schematic  diagram 

of  Vickers  indentation  on  tensile  surface,    (b)  An  elliptical  flaw  on  the  fracture  cross 
section. 
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Furthermore,  we  can  compare  the  obtained  results  from  fracture  surface  analysis  with 
those  from  the  strength-indentation  technique  by  using: 

K,  =  Q.59(E/}if\a,?'"f'.  (4-3) 

Here,  the  fractore  toughness  is  given  in  terms  of  the  critical  intensity  factor, 
Kc,  instead  of  Kj^.  ^ic  is  the  critical  stress  intensity  factor  for  mode  I.  Kjc  is  usually 
determined  from  a  prescribed  fracture  mechanics  test  with  no  residual  stress  associated 
with  the  crack.  K^  is  the  resistance  to  crack  propagation  in  the  presence  of  a  local 
residual  stress.  Kjc  and  Kc  are  expected  to  be  close  in  value.  During  the  test,  the 
mode  one  condition  is  not  guaranteed  because  the  fracture  paths  tend  to  go  to  the 
easiest  cleavage  plane.  Also,  the  stress  condition  on  the  fracmre  surface  is  not  in  a  full 
tension  mode.  Indeed,  sometimes  it  is  a  combination  of  mode  I  and  mode  II.  Thus,  Kc 
is  an  approximation  of  the  mode  I  fracmre  toughness  of  the  material. 

Fractal  Dimension  Determination 

Single  crystal  silicon  fracture  samples  were  carefiiUy  cleaned  and  coated  with 
nickel.  The  samples  were  then  potted  in  epoxy.  In  general,  there  are  two  ways  to  get 
contour  lines:  one  is  to  polish  the  sample  approximately  parallel  to  the  fracmre  surface 
and  the  other  is  to  polish  the  sample  approximately  vertical  to  the  fracmre  surface. 
Thus,  the  so-called  horizontal  contour  (fractal)  dimension  or  vertical  profile  dimension 
can  be  obtained.  The  samples  in  this  smdy  were  polished  approximately  parallel  to  the 
fracmre  surface,  as  shown  in  Fig.  4.3.  As  the  fracmre  surface  is  first  encountered,  a 
section  of  the  fracmre  surface  appears  in  the  polishing  plane.  These  sections  appear  as 
islands  in  the  polishing  plane,  and  are  called  slit-islands.  The  nickel  coating  performs 
two  functions:  It  provides  good  contrast  during  polishing  and  it  helps  to  hold  the 
fracmre  surface  together  during  polishing.    As  polishing  proceeds,  these  islands  begin 
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Fig  4  3  Fractured  samples  are  encapsulated  in  epoxy  and  polished  approximately 

parallel  to  the  fracture  plane.  (Top  View)  Islands  emerge  in  the  polished  plane. 
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to  grow.    The  perimeter  of  the  islands  presents  a  line  or  section  of  a  line  that  can  be 
measured  according  to  Richardson's  equation: 

L(S)  =  k  S^-D,  (4.4) 

where  L(S)  is  the  length  of  a  section  of  a  line  along  the  perimeter  of  the  slit 

island  and  its  value  depends  on  the  measuring  scale,  S; 

S  is  the  measuring  scale  used  to  measure  the  perimeter  and  its  value 

ranged  from  5  to  100  |Lim; 

D  is  the  horizontal  contour  (fractal)  dimension;  and 

k  is  a  constant. 
Richardson  plots  detail  the  change  in  measured  length  of  a  line  as  a  function  of 
scale.    Construction  of  a  Richardson  plot  for  a  fracmre  surface  requires  a  luie  that  is 
representative  of  the  fracture  surface  and  a  range  of  scales  for  measuring  that  line. 

At  the  first  emergence  of  an  island,  polishmg  must  be  carefully  performed.  The 
surface  was  polished  to  a  1  jxm  finish.  Since  the  epoxy  is  transparent,  the  exact 
location  on  the  fracture  surface  can  be  observed,  so  that  measurements  can  be  taken  for 
selected  segments  of  the  perimeter,  e.g.,  in  the  crack  branching  region.  Polaroid 
photographs  were  taken  at  a  magnification  of  400x  and  combined  in  a  montage.  This 
montage  was  then  measured  with  dividers  set  to  various  openings,  as  depicted  in  Fig. 
4.4.  The  length  of  the  section  of  perimeter  was  measured  for  each  divider  setting.  In 
this  way,  the  line  or  perimeter  length  was  computed  as  a  function  of  scales.  A  log-log 
plot  of  length  vs.  scale  results  in  a  line  with  slope  equal  to  1-D. 

The  scale  is  a  measure  of  discernibility.  As  the  scale  becomes  finer  and  finer, 
we  observe  greater  and  greater  detail.  If  the  scale  is  larger  than  the  largest  features, 
the  observation  is  insensitive  to  those  features.  A  curve  will  begin  to  look  fractal  only 
after  the  scale  becomes  smaller  than  such  feamres.  The  scale  lengths  used  in  this  study 
ranged  from  1  unit  length  to  16  unit  lengths  ( 1  unit  length  is  about  5  \xm). 


60 


Fig.  4.4  A  montage  of  the  perimeter  is  measured  with  dividers. 
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M.D.  Simulation  Procedure 

Before  performing  the  MD  simulation,  parameters'^  for  the  potential  and  initial 
conditions  should  be  given  beforehand.  Two  data  files  which  contain  these  data  are 
needed  in  order  to  perform  the  MD  simulation.  One  is  called  the  input  file,  which 
includes  the  initial  conditions  and  the  parameters  of  the  potentials  used.  The  other  one 
is  the  atom-position  file  which  contains  the  initial  position  of  each  Si  atom. 

Several  tasks  were  set  to  be  accomplished  using  MD  simulation.  One  was  to 
compare  the  toughness  of  Si  for  different  orientations.  The  difference  of  stress-strain 
curves  for  different  strain  rates  could  be  obtained  to  compare  the  toughness  as  a 
function  of  strain  rate.  Another  was  to  obtain  the  fractal  dimensions  for  different 
fracture  planes  in  order  to  compare  them  with  the  experimental  results.  Another  goal 
was  to  obtain  the  fracture  strengths  for  each  simulation.  The  comparison  of  fracture 
strength  due  to  different  strain  rates  has  been  performed.  The  comparison  of  fracture 
strength  due  to  different  crack  sizes  has  also  been  performed. 

Determination  of  Input  Data 

Before  performing  the  MD  simulation,  several  conditions  should  be  decided 
first.  All  those  conditions  should  be  illustrated  in  the  input  file  for  the  MD  simulation. 
Those  conditions  are  the  desired  orientation,  temperamre,  strain  rate,  length  of  each 
time  step,  length  of  the  simulation,  and  the  adiabatic  or  isothermal  state. 

At  first,  after  choosing  the  desired  orientation,  a  data  file  which  contains  the 
number  and  the  positions  of  those  atoms  should  be  constructed.  The  data  file  will 
determine  the  initial  position  of  each  atom  for  the  MD  simulation.  In  this  study,  the 
atom  position  for  the  ideal  crystal  structure  will  be  given  at  first  and  the  crystal  will  be 
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allowed  to  thermally  vibrate  for  one  unit  time,  1  pica-second,  before  performing  strain 
pulling. 

The  temperature  to  perform  the  test  should  also  be  given.  The  chosen 
temperature  is  important  because  it  determines  the  energy  state  of  the  system.  Here  the 
temperatore  or  initial  temperature  is  chosen  to  be  300  °K. 

The  strain-rate  should  be  given,  also.  Different  strain-rates  will  give  different 
results.  Basically,  two  categories  will  be  given.  One  is  higher  than  the  speed  of  sound 
while  the  other  is  slower  than  the  speed  of  sound.  Generally,  results  from  different 
strain-rate  experiments  will  be  compared  for  the  same  orientation.  While  for  different 
orientations,  a  strain  rate  of  0.2  was  usually  chosen.  The  reason  for  selecting  0.2  will 
be  explained  in  the  results  and  discussion  section. 

A  careful  choice  of  the  length  of  each  time  step  is  very  important.  It  plays  a 
very  important  role  not  only  on  the  stability  of  the  system  but  also  on  the  length  of  the 
simulation  time.  If  the  length  of  the  time  step  is  chosen  to  be  too  long,  atoms  will 
move  too  far  in  each  time  step.  Thus  the  simulation  will  either  be  less  realistic  or 
become  unstable.  If  the  time  step  is  too  small,  atoms  just  move  a  very  small  distance  at 
each  time  step  and  waste  a  lot  of  computation  time.  For  silicon,  0.005  pica  seconds  is 
a  good  choice  for  the  Stillinger-Weber  and  (ionic)  Coulombic  potentials  because  the 
thermal  vibration  period  is  found  to  be  0.0766  ps  for  Si.^"^  The  length  of  the 
simulation  is  determined  by  the  fracture  of  the  specimen.  For  a  lower  strain-rate  test,  a 
longer  time  should  be  given. 

The  choice  of  an  adiabatic  or  isothermal  condition  is  very  important  to  do  MD 
simulations.  If  the  adiabatic  condition  is  used,  the  temperature  will  increase  several 
thousands  Kelvin  degrees  and  the  system  becomes  unstable  after  fracture.  If  only  the 
Stillinger-Weber  potential  is  used,  the  temperature  after  fracture  is  about  7,000  °K.  If 
the  ionic  Coulombic  potential  is  added  in,  the  temperature  after  fracmre  will  be  about 
50,000°K.   The  reason  is  that  the  ionic  Coulombic  potential  is  always  repulsive.   Thus, 
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the  energy  state  of  the  system  will  increase  as  more  positive  charged  Si  atoms  are 
produced  during  the  fracture  process.  On  the  other  hand,  the  isothermal  condition  will 
result  in  a  more  stable  thermal  state  for  this  study  since  the  temperature  remains 
constant.  Thus,  the  isothermal  condition  was  used  in  this  study.  The  temperature  vs. 
strain  curves  for  both  isothermal  and  adiabatic  condition  are  shown  in  Fig.  4.5  and  4.6. 
The  combinations  of  the  previous  conditions  generally  should  satisfy  different 
needs.  After  these  conditions  have  been  decided,  one  can  begin  to  perform  the  MD 
simulation. 

Simulation  Procedure 

At  first,  toughnesses  for  different  orientations  were  compared.  Three  different 
orientations  were  chosen  in  this  study.  The  fracture  planes  were  chosen  to  be  {100}, 
{110}  and  {111}  as  shown  in  Fig.  4.7.  A  Periodic  boundary  condition  was  applied 
during  the  simulation.  A  constant  expansion  rate  in  the  x-direction  was  applied  for 
each  simulation. 

The  stress-strain  curves  for  different  strain  rates  were  compared  in  order  to  find 
the  dependence  of  strength  on  strain  rate.  The  strain  rates  were  chosen  to  be  0.1,  0.2, 
0.5,  1,  2,  5  x-length/ps. 

In  the  Griffith  criterion  the  strength  of  brittle  materials  depends  on  the  crack 
size  at  which  fracture  begins.  The  introduction  of  a  crack  in  the  ideal  crystal  will  cause 
a  decrease  in  the  fracture  strength.  Larger  cracks  are  expected  to  result  in  lower 
strengths.  A  crack  can  be  introduced  in  two  ways.  One  type  of  crack  simulates  a  void 
by  removing  a  cluster  of  atoms.  The  other  type  simulates  a  planar  sharp  crack  by 
removing  a  layer  of  atoms. 
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Fig.  4.5  The   temperature   vs.    strain   curve   for   isothermal   condition.       The 

temperature  remains  at  about  300°K  before  the  moment  when  fracture  takes  place  at 
strain  »  0.45.   After  then  the  temperature  goes  back  to  about  300  °K. 
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Fig.  4.6  The   temperature    vs.    strain   curve    for    adiabatic    condition.       The 

temperature  remains  at  about  300°K  before  the  moment  when  fracture  takes  place  at 
strain  «  0.45.  After  fracture,  the  temperature  of  the  system  rise  very  fast  up  to  7000° 
K,  then  slows  down. 
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Fig.  4.7  The  {100},  {110},  and  {111}  planes  are  those  fracture  planes  chosen  in 

the  MD  simulation. 
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Fractal  Analysis  Using  Simulation  Results 

After  the  fracture  has  occurred  in  the  simulation  as  shown  in  Fig.  4.8,  an 
important  objective  is  to  obtain  the  generated  fracture  surface  from  the  MD 
simulation. Each  atom  is  assumed  to  be  a  sphere  with  an  electron  cloud  of  a  prescribed 
radius.  The  radius  of  each  atom  is  assumed  to  be  the  radius  of  the  potential  field. 
Thus  the  fracture  surface  appears  as  a  surface  with  intersecting  spheres,  as  shown  in 
Fig.  4.9.  After  obtaining  the  surface,  a  contour  plane  is  used  to  intersect  the  spheres  to 
obtain  a  "slit-island."  The  obtained  slit-island  appears  as  a  plane  with  lots  of  circular 
disks  on  it  as  shown  in  Fig.  4.10.  The  slit-island  technique  as  described  in  the  fractal 
dimension  determination  section  in  the  experimental  procedure  was  applied  here  to 
obtain  the  fractal  dimension.  The  perimeter  of  the  selected  island  was  measured  using 
different  scales  as  described  before  and  the  fractal  dimension  can  be  obtained  from  the 
Richardson  plot. 
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Fig.  4.8  The  stereo-pair  picture  shows  a  frozen  frame  during  fracture. 
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Fig.  4.9  If  each  atom  is  assumed  to  be  as  a  potential  field  with  sphere  shape.   The 

fracmre  surface  will  look  like  a  surface  with  intersecting  spheres. 
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Fig.  4.10         The  slit-island  obtained  from  the  simulated  fracture  surface  looks  like  a 
plane  with  lots  of  circular  disks  on  it. 


CHAPTER  5 
RESULTS  AND  DISCUSSION 


Fracture  toughnesses  and  fractal  dimensions  were  obtained  in  two  ways,  one 
was  from  experimental  measurement  and  another  was  from  molecular  dynamics  (MD) 
simulation.  Fracture  toughnesses  and  fractal  dimensions  were  investigated  for  three 
different  fracture  planes.  This  chapter  discusses  the  toughness  results,  the  fractal 
dimension  results,  and  how  the  fractal  dimensions  for  each  orientation  are  related  to  the 
toughness  results.  This  chapter  is  divided  into  three  sections:  experimental  results, 
sunulation  results  and  comparison  between  measured  and  simulation  results. 
Relationships  between  variables  will  be  discussed  in  each  section. 

Experimental  Results 

In  this  section,  results  obtained  from  the  experimental  measurements  will  be 
given.  These  include  fracture  surface  analysis,  fracture  toughness  measurements, 
fractal  dimension  measurements,  and  the  relationship  between  fractal  dimension  and 
fracture  toughness. 

Fracture  Surface  Analysis 

After  the  specimens  were  broken,  the  fracmre  surfaces  were  examined.  The 
critical  flaw  size  was  measured  to  calculate  the  critical  stress  intensity  factor  using  Eq. 
(4.2).  As  in  Fig.  4.1,  three  fracmre  planes,  i.e.,  {100},  {110}  and  {111},  were  chosen 
and  two  different  tensile  surfaces  in  the  {110}  and  {111}  fracture  planes  were  tested. 


71 


72 


Thus,  there  were  a  total  of  5  different  groups  of  data.   The  fracture  surfaces  are  shown 
in  Figs.  5.1,5.2,5.3.   As  shown  in  Fig.  5.1,  the  fracture  surface  of  the  {100}  fracture 
plane  appears  to  be  the  most  tortuous  one.    It  has  the  smallest  mirror  region  when 
compared  with  the  other  two  orientations.    The  fracture  surface  of  the  {110}  fracture 
plane  with  a  {100}  tensile  surface  has  a  "Batman"-shaped  mirror  region  as  shown  in 
Fig.  5.2(a).     The  frachire  surface  of  the  {110}  frachire  plane  with  a  {110}  tensile 
surface  has  an  inverted  volcano-shaped  mirror  region  as  shown  in  Fig.  5.2(b),  while 
the  frachire  surface  of  the  {111}  fracture  plane  as  shown  in  Fig.  5.3  is  relatively 
smooth  compared  with  the  other  orientations.     "River-marks",  i.e.,  twist  hackle  or 
cleavage  marks,  appear  on  most  of  the  fracture  surfaces  of  this  fracture  plane.    "River- 
marks"  are  twist  hackle  fracture  features  which  appear  to  spread  apart  as  the  crack  path 
moves  away  from  the  origin  of  fracture,  thus  "pointing"  back  to  the  origin.   The  {111} 
frachire  plane  is  the  easiest  cleavage  plane,  as  reported  before. ^^  Shown  in  Fig.  5.3(a) 
is  the  fracture  surface  of  the  {111}  fracmre  plane  with  a  {110}  tensile  surface  while 
that  in  Fig.  5.3(b)  is  the  fracture  surface  of  the  same  frachire  plane  but  with  a  {112} 
tensile  surface.   More  data  are  listed  in  Tables  5.1,  5.2,  and  5.3. 

For  fracture  on  the  same  plane,  there  should  be  no  structural  difference  during 
the  fracmre  process.  However,  the  loading,  i.e.,  strain  would  be  expected  to  be 
different  in  different  locations  due  to  elastic  anisotropy  and  a  changing  strain  field, 
which  m  turn  is  due  to  the  type  of  strain  applied  here,  i.e.,  benduig.  The  fracture 
surface  can  appear  different  due  to  a  change  in  the  tensile  surface  in  the  bending  test. 
Elastic  anisotropy  is  the  reason  the  frachire  surfaces  look  different  in  the  same  fracmre 
plane  without  a  strucmral  change.  Thus  the  fracmre  surfaces  of  the  {110}  fracture 
planes  as  shown  in  Figs.  5.2(a)  and  (b)  look  different  because  the  tensile  surfaces  are 
different.  An  interesting  result  occurs  if  we  rotate  the  fracmre  surface  by  90  degrees 
with  a  {110}  tensile  surface,  as  shown  in  Fig.  5.4;  a  half  mirror  boundary  of  this 
fracmre  surface  is  similar  to  a  half  mirror  boundary  of  the  fracmre  surface  on  the 
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Fig.  5.1  A  typical  fracture  surface  on  the  {100}  fracture  plane. 
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Fig.  5.2  (a)  A  typical  fracture  surface  of  the  {100}  tensile  surface  on  the  {110} 

fracture  plane.  A  "Batman" -like  mirror  is  obvious  on  the  fracture  surface.  And  a 
typical  flaw  is  easily  seen  on  the  fracture  surface,  (b)  A  typical  fracture  surface  of  the 
{110}  tensile  surface  on  the  {110}  fracture  plane.  An  inverse  volcano-like  mirror  is 
obvious  on  the  fracture  surface,  and  a  typical  flaw  is  easily  seen  on  the  fracture 
surface. 
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(a) 


{110} 


(b) 


{112} 


Fig.  5.3  (a)  A  typical  fracture  surface  of  the  {110}  tensile  surface  on  the  {111} 

fracture  plane,  (b)  A  typical  fracture  surface  of  the  {112}  tensile  surface  on  the  {111} 
fracture  plane.  It  is  found  that  both  fracture  surfaces  look  similar  even  though  the 
tensile  surfaces  are  different. 
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Table  5.1     Data  for  Si  samples  fractured  in  the  {100}  plane 


{100}  fracture  plane  ,<11Q>  tensile  surface 

mdeni 

h 

flaw 

load 

Klc* 

KIc** 

load 

width 

thickness 

b 

a 

size 

P 

strength 

0.59(E/ 

1.65* 

No. 

H 

mm 

mm 

mm 

mm 

mm 

lb 

MPa 

MPam^l/2 

MPam^l/2 

al 

0.5 

5.70 

5.35 

a2 

0.5 

6.75 

5.00 

0.068 

0.068 

0.068 

96 

102.5 

1.31 

1.40 

a3 

0.5 

5.85 

4.95 

0.080 

0.080 

0.080 

65 

81.7 

1.10 

1.21 

a4 

0.5 

6.70 

5.40 

0.056 

0.052 

0.054 

118 

108.8 

1.37 

1.23 

average 

1.26 

1.28 

dev 

0.11 

0.08 

bl 

1 

6.25 

5.05 

0.127 

0.094 

0.110 

63 

71.2 

1.18 

1.23 

b2 
b3 

1 

5.90 

5.25 

0.080 

0.080 

0.080 

78 

86.4 

1.37 

1.28 

1 

5.80 

5.05 

0.124 

0.124 

0.124 

54 

65.8 

1.11 

1.21 

b4 

1 

5.90 

5.00 

0.114 

0.072 

0.091 

63 

77.0 

1.25 

1.21 

b5 

1 

5.80 

5.10 

0.127 

0.099 

0.112 

60 

71.7 

1.19 

1.25 

average 

1.23 

1.24 

dev 

0.09 

0.03 

cl 

2 

6.25 

5.20 

0.201 

0.132 

0.163 

55 

58.6 

1.22 

1.23 

c2 

2 

5.90 

5.15 

0.159 

0.131 

0.144 

57 

65.6 

1.32 

1.30 

average 

1.27 

1.27 

dev 

0.05 

0.03 

total 

average 

1.24 

1.26 

total 

dev 

0.09 

0.06 

three-point  bending   span=27mm 

*  KIc  is  calculated  from  Eq.  (4.3) 

**KIc  is  calculated  from  Eq.  (4.2) 
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Table  5.2     Data  for  Si  samples  fractured  in  the  {1 10}  plane 


{100}  tensile  surface 

indenl 

t 

flaw 

load 

Klc* 

Klc** 

No. 

load 

width 

thickness 

b 

a 

size 

P 

strength 

0.59(E/ 

1.65* 

l^g 

mm 

mm 

um 

um 

um 

N 

MPa 

MPam'^1/2 

MPam'^1/2 

1 

0.9 

7.44 

7.03 

66.2 

84.8 

74.0 

322.6 

103.1 

1.46 

1.46 

2 

0.9 

7.31 

7.02 

70.2 

6^.3 

68.0 

0.0 

107.4 

1.51 

1.46 

i 

0.9 

7.31 

7.00 

66.2 

79.5 

72.0 

249.2 

81.8 

1.23 

1.14 

4 

0.9 

7.31 

6.99 

79.5 

92.8 

86.0 

262.5 

86.7 

1.29 

1.32 

:> 

0.9 

7.53 

7.15 

95.4 

98.0 

96.5 

266.9 

81.5 

1.23 

1.31 

b 

0.9 

6.87 

7.16 

92.8 

90.0 

91.5 

238.0 

79.4 

1.21 

1.25 

/ 

0.9 

6.58 

7.03 

94.0 

94.0 

195.8 

70.8 

1.10 

1.13 

8 

0.9 

6.89 

7.01 

97.0 

97.0 

242.6 

84.2 

1.26 

1.35 

y 

0.9 

7.07 

7.07 

102.0 

102.0 

202.4 

67.3 

1.06 

1.12 

10 

0.9 

6.81 

7.39 

95.0 

89.8 

89.8 

233.6 

73.9 

1.14 

1.15 

11 

12 

0.9 

7.55 

6.95 

89.8 

92.4 

198.0 

63.9 

1.02 

1.01 

1.4 

7.11 

7.37 

130.0 

143.0 

136.3 

164.6 

50.1 

0.95 

0.95 

u 

1.4 

6.88 

7.11 

136.0 

0.87 

14 

1.4 

7.84 

7.16 

132.5 

124.0 

127.9 

193.5 

56.6 

1.05 

1.05 

!:> 

1.4 

7.07 

7.44 

119.0 

106.0 

112.1 

229.1 

68.8 

1.21 

1.23 

lb 

1.8 

6.83 

6.92 

126.0 

172.0 

147.0 

198.0 

71.1 

1.32 

1.41 

1/ 

1.8 

7.43 

7.02 

242.5 

77.8 

1.41 

18 

1.8 

7.82 

6.94 

220.0 

220.2 

68.7 

1.29 

19 

1.8 

6.90 

7.35 

185.0 

185.0 

200.2 

63.2 

1.21 

1.41 

20 

0.5 

6.87 

7.28 

63.4 

52.8 

57.8 

269.2 

87.0 

1.12 

1.09 

21 

0.5 

7.36 

6.60 

52.8 

52.8 

52.8 

262.5 

96.3 

1.20 

1.15 

22 

2 

7.57 

6.90 

132.0 

142.6 

137.2 

204.3 

66.8 

1.29 

1.29 

2i 

2 

7.70 

7.32 

190.0 

158.4 

173.5 

191.3 

54.5 

1.11 

1.18 

24 

3 

7.41 

7.24 

195.4 

142.6 

166.9 

209.1 

63.3 

1.37 

1.35 

2^ 

3 

6.87 

6.95 

169.0 

184.8 

176.7 

186.9 

66.2 

1.42 

1.45 

2b 

5 

6.29 

7.15 

322.1 

294.6 

307.6 

133.5 

48.8 

1.28 

1.41 

11 

5 

7.45 

7.37 

316.8 

253.4 

283.3 

169.1 

49.1 

1.29 

1.36 

avg. 

1.19 

1.23 

dev. 

0.08 

0.08 

{110}  tensi 

le  surface 

indent 

flaw 

load 

Klc* 

Klc** 

No. 

load 

width 

thickness 

b 

a 

size 

P 

strength 

0.59(E/ 

1.65* 

Kg 

mm 

mm 

um 

um 

um 

N 

MPa 

MPam'^1/2 

MPam'^1/2 

41 

0.9 

7.35 

6.97 

158.4 

105.6 

129.3 

160.2 

52.76 

0.89 

0.98 

42 

0.9 

7.55 

7.39 

95.0 

95.0 

95.0 

278.1 

79.30 

1.20 

1.26 

4i 

0.9 

7.10 

7.56 

89.8 

100.3 

94.9 

233.6 

67.71 

1.07 

1.07 

44 

0.9 

6.93 

7.84 

105.6 

105.6 

105.6 

229.1 

63.24 

1.01 

1.06 

4.b 

0.9 

6.76 

7.03 

84.5 

79.2 

81.8 

204.7 

72.06 

1.12 

1.08 

46 

0.9 

7.44 

7.26 

84.5 

95.0 

89.6 

215.8 

64.73 

1.04 

1.01 

4y 

5.0 

7.08 

7.21 

380.1 

337.9 

358.4 

100.1 

31.99 

0.93 

0.99 

48 

5.0 

7.39 

7.26 

269.3 

316.8 

292.1 

142.4 

43.00 

1.17 

1.21 

avg. 

1.05 

1.07 

dev. 

0.10 

0.09 

*  Klc  is  C3 

ilculate< 

1  from  the 

liq.  (4 

.3) 

**  Klc  is  cal( 

:ulated  fn 

Dm  the  Eq.(4 

.2) 
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Table  5.3     Data  for  Si  samples  fractured  in  the  {111}  plane 


110}  tensile  surface 

indent 

flaw 

loac 

KIc* 

KIc** 

load 

width 

1  thickness 

b 

a 

size 

P 

strength 

0.59(E/ 

1.65* 

No. 

J^§ 

mm 

mm 

mm 

mm 

mm 

lb 

MPa 

MPam'^1/2 

MPam^l/2 

cl 

6.15 

4.00 

0.082 

0.092 

0.087 

33 

76.1 

1.20 

1.17 

c2 

5.80 

4.95 

0.093 

0.102 

0.098 

46 

73.4 

1.17 

1.20 

c3 
Al 

6.30 

5.05 

0.091 

0.103 

0.097 

40 

56.5 

0.96 

0.92 

6.15 

5.60 

0.110 

0.112 

0.111 

82 

70.9 

1.14 

1.23 

A2 

6.20 

5.20 

0.103 

0.124 

0.113 

78 

77.6 

1.22 

1.36 

avg 

1.14 

1.18 

dev 

0.09 

0.15 

dl 

0.7 

6.15 

4.65 

0.078 

0.074 

0.076 

44 

75.1 

1.09 

1.08 

d2 

0.7 

6.20 

5.15 

0.058 

0.063 

0.061 

72 

99.4 

1.34 

1.28 

d3 

0.7 

6.10 

4.90 

0.080 

0.086 

0.083 

48 

74.4 

1.08 

1.12 

d4 

0.7 

6.25 

5.05 

0.065 

0.085 

0.075 

56 

79.7 

1.14 

1.14 

avg 

1.16 

1.15 

dev 

0.11 

0.07 

el 

2 

6.30 

4.80 

0.171 

0.200 

0.185 

34 

53.2 

1.09 

1.19 

e2 

2 

6.25 

4.75 

0.149 

0.155 

0.152 

40 

64.4 

1.26 

1.31 

e3 

2 

5.70 

4.70 

0.172 

0.173 

0.172 

29 

52.3 

1.08 

1.13 

e4 

2 

5.95 

4.95 

0.160 

0.118 

0.137 

38 

59.1 

1.18 

1.14 

avg 

1.15 

1.19 

dev 

0.07 

0.07 

{112}tensi] 

e  surface 

fl 

2 

6.15 

6.25 

0.108 

0.115 

0.111 

68 

64.2 

1.26 

1.12 

12 

2 

5.35 

6.20 

0.125 

0.144 

0.135 

60 

66.2 

1.29 

1.27 

13 

2 

5.25 

6.15 

0.165 

0.125 

0.143 

50 

57.1 

1.15 

1.13 

f4 

2 

5.80 

6.10 

0.165 

0.183 

0.174 

48 

50.5 

1.05 

1.10 

avg 

1.19 

1.15 

dev 

0.09 

0.07 

§2 

6.40 

5.80 

0.079 

0.094 

0.086 

67 

70.6 

1.14 

1.08 

5.50 

6.20 

0.105 

0.105 

0.105 

70 

75.1 

1.19 

1.27 

g3 

5.85 

5.65 

0.086 

0.099 

0.092 

52 

63.2 

1.05 

1.00 

Bl 

6.30 

6.05 

0.109 

0.098 

0.103 

101 

73.1 

1.17 

1.23 

B2 

5.80 

5.70 

0.086 

0.071 

0.078 

93 

82.3 

1.28 

1.20 

avg 

1.16 

1.16 

dev 

0.07           0.10      1 

*  KIc  is  calculated  from  Eq.  (4. 

3) 

**KIc  is  calculated  from  Eq.  (4.2) 

A,B    span  =  25.5  mm 

c,d,e,f,  and  g 

span  =34 

three-point  bending 
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(a) 


{110} 


(b) 


{100} 


Fig.  5.4  An  interesting  result  is  that  if  we  rotate  the  frachire  surface  of  the  one 

with  a  {110}  tensile  surface  90  degrees,  as  shown  in  (a),  then  a  half  mirror  boundary  of 
this  fracture  surface  is  similar  to  half  a  mirror  boundary  of  the  fracture  surface  on  the 
{110}  fracture  plane  with  a  {100}  tensile  surface  without  rotation,  as  shown  in  (b). 
This  implies  that  a  fracture  mirror  boundary  for  a  specimen  fractured  in  tension  from 
an  internal  defect  would  mimic  the  trace  of  the  elastic  constant  values  in  that  plane. 
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{110}  fracture  plane  with  a  {100}  tensile  surface  without  rotation.^^  xj^g  implies  that 
a  fracture  mirror  boundary  for  a  specimen  fractured  in  tension  from  an  internal  defect 
would  mimic  the  trace  of  the  elastic  constant  values  in  that  plane. 

Freiman  et  al.52  found  that  the  fracture  mirror  shapes  in  sapphire  for  pseudo- 
cleavage  type  fracture  have  much  more  complex  geometry,  as  shown  Fig.  5.5,  than  the 
similar  fractures  in  glasses  which  yield  more-or-less  circular  fracture  mirrors.  They 
relate  the  mirror  boundaries  to  the  critical  stress  intensity  factor.  The  complex  mirror 
shapes  investigated  in  sapphire  implies  that  the  critical  stress  intensity  factor  along  the 
different  crystal! ographic  directions  of  the  mirror  boundary  is  not  constant.  Young's 
modulus  also  plays  an  important  role  for  the  complex  shape  of  mirror  boundary.  The 
Young's  modulus  normal  to  the  crack  plane  is  a  constant.  However,  the  Young's 
moduli  either  parallel  or  perpendicular  to  the  mirror  boundary  in  the  crack  plane  can 
vary  as  one  traverses  around  the  boundary. ^^  As  seen  in  Fig.  5.5,  The  Young's 
moduli  for  directions  parallel  and  perpendicular  to  the  crack  front/mirror  boundary 
plotted  as  a  function  of  angle  of  rotation  about  the  tensile  axis  possesses  maxima- 
minima  patterns  which  are  similar  to  the  fracture  mirror  geometry.  A  similar  result  for 
single  crystal  silicon  is  expected  and  a  schematic  of  the  expected  mirror  boundary  is 
plotted  in  Fig.  5.5. 

In  contrast  to  the  fracture  surface  of  the  {110}  fracture  plane,  the  fracture 
surfaces  of  the  {111}  fracture  plane  look  the  same  even  when  the  chosen  tensile 
surfaces  are  different  as  shown  in  Fig.  5.3.  Most  probably,  the  fracture  surface  of  the 
easiest  cleavage  plane  will  not  be  affected  by  the  chosen  tensile  surface  because  the 
crack  growth  will  remain  on  the  easiest  cleavage  fracture  plane  as  the  lowest  energy 
path.  Also,  the  size  of  the  specimen  is  not  large  enough  to  show  crack  branching  at 
these  tested  fracture  stress  levels.  Thus,  the  fractography  would  appear  the  same  no 
matter  what  the  change  in  the  tensile  surface  at  these  test  specimen  sizes.  The  available 
fracture  paths  for  the  specific  fracture  planes  are  another  important  factor  to  determine 
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(a) 


t  <i2io> 


(b) 


Fig.  5.5  (a)  Anisotropy  mirror  shapes  on  fracture  surfaces   in  single  crystal 

alumina;  the  corresponding  relative  magnitudes  of  Young's  modulus  for  directions  ( in 
the  plane)  parallel  and  perpendicular  to  the  local  mirror  boundary  are  given.  Larger 
values  of  E  are  reflected  as  greater  distances  from  the  origin,  (b)  A  schematic  of  the 
expected  mirror  boundary  for  the  {110}  fracture  plane  is  plotted. 
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the  topography  of  the  fracture  surface.    A  more  detailed  discussion  of  macroscopic 
crack  branching  follows. 

The  crack  branching  plane  is  an  interesting  topic  in  this  study.  From  the 
observation  of  the  obtained  fracture  surfaces,  the  fracture  surface  of  the  {111}  fracture 
plane  is  relatively  smooth  such  that  no  crack  branching  is  observed  as  shown  in  Fig. 
5.6.  For  fracmre  in  the  {110}  plane,  the  crack  either  doesn't  branch  or  tends  to  branch 
to  the  available  {111}  planes  as  shown  in  Figs.  5.7  and  5.8.  In  Fig.  5.7,  the  crack 
growth  (shown  by  shaded  regions)  occurs  on  the  {110}  plane  first  and  deviates  to 
{111}  plane  as  shown  in  (c)  and  (d)  of  Fig.  5.7.  In  Fig.  5.8,  the  crack  growth  (shown 
by  shaded  regions)  occurs  on  the  {110}  plane  first  and  deviates  to  {111}  planes  as 
shown  in  (c)  and  (d)  of  Fig.  5.8. 

Crack  branching  for  the  {100}  fracture  plane  is  complex.  As  seen  in  Figs.  5.9 
and  5.10,  the  system  of  crack  branching  angles  varies  for  one  fracture  plane  as  well  as 
for  the  loading  direction  changes,  e.g.  from  the  {100}  tensile  surface  to  the  {110} 
tensile  surface. 

For  those  fracture  surfaces  of  the  {100}  fracture  plane  with  a  {100}  tensile 
surface,  the  crack  easily  branches  to  the  {111}  plane  as  shown  in  Fig.  5.9.  If  the 
tensile  surface  is  chosen  to  be  {110},  the  crack  will  remain  on  the  {100}  plane  for 
about  2  to  3  mm  before  branching  to  the  {111}  plane  if  the  indentation  load  is  larger 
than  2  Kgw,  as  shown  in  Fig.  5.10(b).  If  the  indentation  load  is  smaller  than  2  Kgw, 
the  crack  will  branch  in  a  distance  smaller  than  2mm  as  shown  in  Fig.  5.10(c). 

It  is  also  found  that  the  fracmre  mirror  region  is  not  really  very  smooth.  ^^  The 
reason  that  we  say  it  is  "mirror-like"  is  because  the  tortuosity  of  the  region  is  smaller 
than  the  wave-length  of  visible  light  as  shown  in  Fig.  5.11.  Indeed,  if  a  scaiming 
tunneling  microscope  (STM)  is  used  to  observe  the  mirror  region,  then  the  surface  is 
not  smooth  at  the  scale  at  which  the  STM  is  used.  12  As  in  Fig.  5.12,  the  surface  on 
the  mirror  region  is  quite  rough  at  that  magnification.     Thus,  the  geometry  in  the 
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^         {110} 


{112} 


(b) 


{110} 


_      {1111 


(c) 


{112} 


{111)       _ 


{110} 
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Fig.  5 .6  Crack  path  for  the  {111}  fracture  plane. 
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(d) 


{100} 


{110} 
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{100} 


{110} 


Fig.  5.7  The  crack  path  for  the  {110}  fracture  plane  with  a  {100}  tensile  surface. 

The  crack  growth  (shown  by  shaded  regions)  occurred  on  the  {110}  plane  first  and 
deviated  to  the  {111}  plane  as  shown  in  (c)  and  (d).  (a)  shows  a  schematic  of  the 
fracture  bar.   (b),  (c),  and  (d)  show  the  different  modes  observed. 
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(a) 


(b) 


^  {110} 

{100} 


{110} 


{100} 


{110} 
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/     \ 

{110}/ 

{100}              .^ 

''{110} 
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(d) 


{111} 


Fig.  5.8  The  crack  path  for  the  {110}  fracture  plane  with  a  {110}  tensile  surface. 

The  crack  growth  (shown  by  shaded  regions)  occurred  on  the  {110}  plane  first  and 
deviated  to  the  {111}  planes  as  shown  in  (c)  and  (d).  (a)  shows  a  schematic  of  the 
fracture  bar.  (b),  (c),  and  (d)  show  the  different  modes  observed. 
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(a) 


-A 


{100} 


{100} 


(b) 


Fig.  5.9  The  crack  paths  (shown  by  shaded  regions)  for  the  {100}  fracture  plane 

with  a  {100}  tensile  surface,     (a)  shows  a  schematic  of  fracture  bar.     (b)  shows 
different  modes  observed. 
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(a) 


■^ 


{110} 


(b) 


(c) 


Fig.  5.10  The  crack  paths  (shown  by  shaded  regions)  for  the  {100}  fracture  plane 
with  a  {110}  tensile  surface,  (a)  shows  a  schematic  of  fracture  bar.  If  the  tensile 
surface  is  chosen  to  be  {110},  the  crack  will  remain  on  the  {100}  plane  for  about  2  to  3 
mm  before  branching  to  the  {111}  plane  if  the  indentation  load  is  larger  than  2  Kgw  as 
shown  in  (b).  If  the  indentation  load  is  smaller  than  2  Kgw,  the  crack  will  branch  in  a 
distance  smaller  than  2  mm  generally  as  shown  in  (c). 
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Fig.  5.11  Schematic  of  vertical  profile  of  the  mirror,  mist  and  hackle  regions  of  a 

fracture  surface.  The  reason  why  the  mirror  region  is  mirror-like  is  because  the  tortuosity 
of  the  region  is  smaller  than  the  wave-length  of  visible  light  as  shown  above. 
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Fig.  5.12         The  fracture  surface  obtained  from  scanning  tunneling  microscopy. 
(Courtesy  of  Mitchell  and  Bonnell^^) 
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mirror  region  appears  to  be  smooth  to  the  unaided  eye,  but  appears  to  be  tortuous  if  the 
STM  is  used  as  the  tool  of  observation. 

The  reason  that  the  fractographic  features  described  above  are  important  is  that 
the  distances  to  crack  branching  have  been  shown  to  be  related  to  different  energy 
levels  at  branching  and  the  fracture  toughness  of  the  material. ^^  In  addition,  fracture 
surface  features  have  also  been  shown  related  to  the  fractal  dimension  of  the  surface. ^-^ 
Thus,  the  fracture  surface  features  are  intimately  coimected  to  the  fracture  process  as 
described  by  the  fracture  toughness,  branching  constants  and  the  fractal  dimensions. 
These  inter-relationships  will  be  discussed  later  in  this  chapter. 

Fracture  Toughness  Measurement 

Single  crystal  Si  is  anisotropic  and  the  fracture  toughness  varies  with  the 
fracture  plane  (and  the  tensile  surface  if  a  bending  technique  is  used  to  fracture  the  Si). 
Chen  and  Leipold^^  determined  the  toughness  for  single  crystal  silicon  tested  in 
different  orientations.  They  found  that  the  {111}  fracture  plane  is  the  easiest  cleavage 
plane  and  the  {100}  fracture  plane  is  the  toughest  one.  They  used  an  indentation 
technique  to  control  the  crack  sizes  on  the  tested  specimens  before  fracture.  However, 
in  their  work,  they  did  not  consider  the  residual  stress  caused  by  the  indentation.  The 
equation  which  they  used  to  calculate  the  critical  stress  intensity  doesn't  consider 
residual  stress: 

Kc  -  aMB(7ra/Q)l/2,  (5.1) 

where  K^  is  the  critical  stress  intensity  factor, 

a  is  the  maximum  tensile  stress, 

a  is  the  flaw  size, 

Mg  is  the  elastic  stress  intensity  magnification  factor,  and 
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Q  is  the  flaw  shape  parameter  calculated  from  the  elliptical  integral  of 
the  second  kind. 

From  the  given  flaw  size  data  in  their  paper,22  the  Mg  value  was  calculated  to 
be  1.03  and  Q  was  approximately  1.47.  After  substituting  these  data  into  equation 
(5.1),  the  equation  can  be  approximately  rewritten  as 

Kc=1.24aVc.  (5.2) 

If  we  modify  Eq.  (5.2)  to  include  the  effect  of  residual  stress  (due  to  the  indentation 
process),  the  new  modified  equation  becomes  Kc  =  1.65afC       i.e.,  Eq.  (4.2).    ' 
Their  data  using  Eq.  (4.2)  fall  in  the  same  range  as  those  measured  in  this  study.  The 
measured  data  and  modified  data  from  Chen  and  Leipold  are  presented  in  Table  5.4. 

Table  5.4  Comparison  of  toughness  values  using  different  techniques 


fracture  plane 

tensile  surface 

S.I. 
toughness* 
(  MPaV^) 

F.S.A. 
toughness** 

(  MPaVm) 

literature 
toughness 
(  MPaV^) 

{100} 

{110} 

1.24  ±0.09 

1.26  ±0.06 

1.26  ±0.07 

{110} 

{100} 
{110} 

1.19  ±0.08 
1.05  +  0.10 

1.23  ±0.08 
1.07  ±0.09 

1.19  ±0.13 

{111} 

{110} 
{112} 

1.19  +  0.10 
1.21  ±0.09 

1.17  ±0.09 
1.16  ±0.08 

1.09  ±0.09 

*Calculated  from  the  strength  indentation  (S.I.)  technique, 

,1/8 


Kc=0.59 


f- 


OfV-) 


,1/3 


3/4 


1/2 


**Calculated  from  fracture  surface  analysis  (F.S.A.),  K^  =  1.65afC 
%y  Chen  and  Leipold. ^^  Originally  these  values  were  calculated  without  considering 
residual  stress  and  were  found  to  be  0.95+0.05  ({100}  plane),  0.90±0.11  ({110} 
plane),  and  0.82±0.07  ({111}  plane)  MPaVm. 
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Gilman^,  Jaccodine,^^  Myers  &  Hillsberry,^^  St.  Johii,^^  and  Kalwani^^  also 
investigated  the  fracture  toughness  in  single  crystal  silicon.  Those  data  and  applied 
techniques  are  listed  in  Appendix  A. 

Single  crystal  Si  is  a  diamond  cubic  structure  which  is  considered  to  be  one  type 
of  face  centered  cubic  structure.  A  simple  approach^^  which  has  been  accepted  to 
estimate  the  surface  energy  of  different  planes  in  a  FCC  structure  is  to  calculate  the 
number  of  dangling  bonds  on  the  specific  plane  for  the  FCC  materials.  The  surface 
energy  can  be  obtained  using  the  equation: 

Y=-s-N{hkl}-n{hkl}.  (5.3) 

where  8  is  the  energy  needed  to  fracture  the  Si-Si  bond, 

—  is  used  because  two  surfaces  are  generated, 

N{hkl}  is  the  number  of  atoms  per  area  on  the  {hkl}  plane,  and 
ii{hkl}  is  the  number  of  dangling  bonds  per  atom  on  the  {hkl}  plane. 
By  using  the  given  equation,  we  can  find  the  surface  energy  for  the  FCC 
structure  on  the  {100},  {110},  and  {111}  planes  separately  as  shown  in  Fig.  5.13. 
Since  the  number  of  atoms  on  the  {100}  plane  is  1  +  1/4  ■  4  as  shown  in  Plot  (b), 
which  is  2,  and  the  area  A  where  those  atoms  lie  is  a  •  a  where  a  is  the  length  of  the 
side  of  a  FCC  unit  cell,  Nngo}  is  found  to  be  lid?-.  Since  the  number  of  dangling 
bonds  per  atom  on  the  {100}  plane  is  4  as  shown  in  Plot  (c),  y(ioo}  is  found  to  be 

12  12 

—  ■£--y4.  In  the  same  way,  Ynio}  is  found  to  be  —•£■  -5,  and  Y{iii}  is  found  to 

■^       a  -^       -SI  ZsL 

be  Y^-J^^-^-    Thus  y|ioo}  :  Y{110}  :  y{lll}  =  8  :  7.07  :  6.92  =  1.16  :  1.02  :1. 

This  ratio  means  that  the  {111}  plane  for  the  FCC  structure  has  the  lowest  surface 
energy  and  the  {100}  plane  has  the  highest  surface  energy. 
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Fig.  5.13         The  surface  energy  for  FCC  structures  on  the  (100),  (110),  and  (111) 
planes  can  be  found  using  a  dangling  bond  calculation. 
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However,  single  crystal  Si  has  a  diamond  cubic  (DC)  structure  which,  in  fact,  is 
different  from  FCC  because  it  is  not  closed  packed.  We  can  find  the  surface  energy  for 
the  DC  structure  on  the  {100},  {110},  and  {111}  planes  separately  as  shown  in  Fig. 
5.14.  However,  the  DC  structure  is  not  a  close-packed  structure,  and  the  dangling 
bonds  will  not  be  linked  to  all  neighbors.  Only  some  specific  neighbors  in  the  first 
shell  will  share  the  dangling  bond.  The  previous  method  cannot  be  completely 
followed  and  some  modifications  must  be  made.  As  shown  in  the  figure,  the  dangling 
bond  is  expressed  as  a  Ime.  The  surface  energy  can  be  obtained  using  the  equation: 

Y=-£-n{hkl}--.  (5.3) 

where  e  is  the  energy  needed  to  fracture  the  Si-Si  bond, 

-  is  used  because  two  surfaces  are  generated, 

2 

nfhki}  is  the  number  of  dangling  bonds  on  the  {hkl}  plane,  and 
A  is  the  area  where  those  atoms  lie  on. 
Since  the  number  of  dangling  bonds  on  the  {100}  plane  is  4  as  shown  in  Plot 
(c),  and  the  area  A  where  those  atoms  lie   is  a  •  a  and  a  is  the  length  of  the  side  of  a 

silicon  DC  unit  cell.   Y{ioo}  is  found  to  be  --s-^.    If  s  =  3.5-10-12  erg/pair-atom22 

and  a  =  5.43  A  is  substimted  into  the  equation,  yjioo}  is  found  to  be  2.3  J/m^.    (The 
Si-Si  bond  energy  is  found  to  be  varied  from  3.2  to  3.6  erg/pair-atom.24,64,68,69)    in 

1  4 

the  same  way,  y|iio}  is  found  to  be  --s---^^  =  1.67  J/m^,  and  yim}  is  found  to  be 

•^       V^a 

1  2 

-s ^ —   =  1.36  J/m2.    These  values  are  close  to  the  experimental  value,  2.1 

2  V3a2/2 

J/m^,  which  is  obtained  using  double  cantilever  beam  tests. ^^'^^    Thus  Y{ioo}  :  Y{iio}  ■ 
y{lll}  =  2.3  J/m2  :  1.67  J/m2  :  1.36  J/m2  =  1.69  :  1.23  :  1.    So,  from  dangling  bond 
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Fig.  5.14         The  surface  energy  for  DC  structures  on  the  (100),  (110),  and  (HI) 
planes  can  be  found  using  a  dangling  bond  calculation. 
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calculations  for  the  DC  structure,  the  {111}  plane  has  the  lowest  surface  energy  and  the 
{100}  plane  has  the  highest  surface  energy,  just  as  before. 

A  direct  measurement  of  surface  energy  of  single  crystal  silicon  using  a 
cleavage  technique  was  done  by  Jaccodine.64  The  measured  surface  energy  depends  on 
determining  the  force  necessary  to  just  move  a  crack  that  is  already  present  in  the 
specimen  under  test.  Determining  this  force  and  the  lever  arm  to  the  crack  tip,  along 
with  other  specimen  geometry  measurements,  allows  one  to  calculate  the  surface 
energy.  Jaccodine  calculated  the  strength  of  the  Si-Si  bond  to  be  45.5  Kcal/mole, 
which  is  equal  to  3.2-10-12  erg/pair-atom.  From  this  value  and  the  bond  density,  the 
surface  energies  for  different  planes  were  found  to  be  2.1  J/m2  for  the  {100}  plane,  1.5 
J/m2  for  the  {110}  plane,  and  1.2  J/m^  for  the  {111}  plane. 

Hesketh  et  al.69  also  calculated  the  surface  free  energy  for  single  crystal  silicon 
in  different  planes.  Using  the  unit  cell  dimension  of  5.432  A  and  a  Si-Si  bond  energy 
of  2.9  10-12  erg/pair-atom,  the  surface  energies  for  different  planes  were  found  to  be 
2.0  J/m2  for  the  {100}  plane,  1.4  J/m2  for  the  {110}  plane,  and  1.2  J/m2  for  the  {111} 
plane. 

A  table  of  comparison  for  the  surface  energy  obtained  by  the  different  authors 
are  listed  in  Table  5.5.  Any  difference  in  values  is  due  to  the  variance  of  Si-Si  bond 
energy  obtained  from  the  different  techniques.  However,  the  values  are  very  similar. 

Thus,  dangling  bond  calculations  show  that  the  {100}  plane  is  the  most  difficult  f 

cleavage  plane  while  the  {111}  plane  is  the  easiest  cleavage  plane.  The  trend  for  the 
toughness  from  dangling  bond  calculations  is  the  same  as  experimental  measurements, 
but  the  numerical  ratio  for  different  fracture  planes  from  the  dangling  bond  calculations 
is  not  the  same  as  the  experimentally  measured  results. 


t 
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Table  5.5  Comparison  of  surface  energy  obtained  from  different  authors. 


fracture  plane 

Si-Si  bond 
(erg/pair-atom) 

{100} 
(J/m^) 

{110} 

(J/m2) 

{111} 

(J/m^) 

this  work 

3.5  •  10-12 

2.3 

1.7 

1.4 

Jaccodine 

3.2  ■  10-12 

2.1 

1.5 

1.2 

Hesketh  et  al.^^ 

2.9  •  10-12 

2.0 

1.4 

1.2 

From  dangling  bond  calculations,  the  frachire  energy  on  the  {100}  plane  is 
almost  twice  as  high  as  that  on  the  {111}  plane.  The  reason  may  be  that  the  dangling 
bond  calculations  do  not  consider  the  angle  between  the  bonds  and  the  specific  fracture 
plane.  For  example,  the  angle  between  the  dangling  bond  and  the  {111}  plane  is  90° 
while  the  angle  between  the  dangling  bond  and  the  {100}  plane  is  about  40°.  If  the 
angle  between  the  bond  and  fracture  plane  is  90°,  more  energy  is  needed  when  two 
atoms  are  pulled  a  distance  dx  away  from  each  other  as  compared  to  the  case  for  an 
angle  less  than  90°.  An  example  is  illustrated  in  Fig  5.15(a),  in  which  two  atoms  are 
pulled  a  distance  dx  in  the  {111}  plane.  The  energy  consumed  for  this  case  is  1/2-K- 
dx'^  since  the  angle  between  the  force  and  displacement  is  0°,  where  K  is  the  force 
(spring)  constant  between  the  two  atoms.  But  for  the  {100}  plane,  as  shown  in  Fig. 
5.15(b),  the  angle  between  the  force  and  the  displacement  is  not  0°.  Thus  if  two  atoms 
were  pulled  a  distance  dx  away  from  each  other,  the  achial  movement  (distance) 

The  quantitative  relation 


between  the  two  atoms  is     ^(x  +  dx)^ +y^ -^x^ -i-y^ 


between  dx  and   ^j{x  +  dx)'^  +y^  -^jx^  +y^    can  be  determined  by  the   following 
procedure: 
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displacement  between  two  atoms  =  -(/(x  +  dx)^  +  y^  -  V  ^^  +  Y^  <  dx 


1  y 

consumed  energy  <  —  ■  E  •  dx 

2 


Fig.  5.15         The  energy  needed  to  elongate  the  Si-Si  bond  is  dependent  on  the  angle 
between  the  dangling  bond  and  the  pulling  direction. 
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V(x  +  dx)2+y2-Vx2+y2         ?        dx.  (5.4) 


If  we  adding  V^^  +  y^  at  both  sides,  then  the  above  equation  changes  to 

^/(x  +  dx)2+y2  ?  ^x2+y2   +dx.  (5.5) 

If  we  square  both  sides,  the  above  equation  changes  to 

x^  +  dx^  +  y^  +  2xdx         l  x^  +  y^  +  dx^  +  2■dx■^J^?+y^ .  (5.6) 

If  we  eliminate  x^  +  dx^  +  y^  on  both  sides,  then  the  l^  mark  can  be  found  to  be 

"  <  " ,  and  the  relation  for  the  above  equation  is 


^(x  +  dxf+y^  -^+y^     <         dx.  (5.7) 

Eq.  (5.7)  means  that  the  energy  needed  to  pull  a  section  dx  in  the  {100}  plane  is 
smaller  than  l/l-K-dx"^,  as  for  the  {111}  plane.  Thus,  if  only  one  bond  is  considered, 
the  energy  consumed  to  pull  one  bond  in  the  {100}  plane  is  smaller  than  that  for  the 
{111}  plane  for  the  same  pulling  distance.  Although  the  ratio  of  toughness  between 
{100}  and  {111}  is  calculated  to  be  2.3/1.36,  the  actual  value  should  be  less  than  that. 

In  the  same  way,  if  dx  is  pulled  in  the  x  direction,  the  elongation  along  the 
dangling  bond  is  longer  for  atoms  laid  on  the  {111}  plane  than  those  laid  in  the  {100} 
plane  since  the  {111}  plane  is  perpendicular  to  the  x  direction.  If  a  bond  will  break 
after  a  specific  elongation,  then  when  one  expands  the  single  crystal  Si  perpendicular  to 
the  {100}  plane  (or  along  the  [100]  direction),  more  deformation  is  needed  to  fracture 
the  specimen. 

Thus,  when  we  observe  the  stress-strain  curve  for  the  {111},  {HO},  {100} 
orientation,  we  can  find  that  the  fractore  strain  for  the  {100}  plane  is  the  highest.  (The 
resuU,  which  will  be  discussed  later,  is  the  same  as  the  stress-stram  curve  obtained 
from  the  MD  simulation,  which  shows  the  fracmre  strain  for  the  {100}  fracture  plane  is 
larger  than  the  other  two  orientations.)  Although  the  modulus  in  the  {100}  plane  is 
lower  than  the  other  two  planes,  the  toughness  (here  the  toughness  is  defined  as  the 
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area  under  the  stress  strain  curve)  for  the  {100}  plane  will  be  higher  than  those  of  the 
other  two  planes  based  on  surface  energy  calculation  and  bond  angle. 

Fractal  Analysis 

The  fractal  dimension  for  each  orientation  is  measured  using  several  different 
slit-island  lines.  Generally,  about  8  or  more  different  sets  of  slit-islands  for  each 
orientation  are  obtained.  A  section  of  the  perimeters  were  photographed  using  an 
optical  microscope  at  400X.  Different  scales  are  used  to  measured  the  length  of  a 
portion  of  the  perimeter  of  the  obtained  slit-island  as  shown  in  Fig.  2.5.  For  each  slit- 
island,  two  sets  of  data  of  length  vs.  scale  are  obtained.  Each  set  can  be  used  to  find 
one  fractal  dimension  using  Richardson's  technique  (Eq.  4.4).  A  log-log  plot  of  length 
vs.  scale  results  in  a  line  with  a  slope  equal  to  1-D  as  shown  in  Fig.  2.6.  The 
measured  fractal  dimension  is  determined  valid  if  the  correlation  coefficient  (r2)  is 
higher  than  0.8.  The  fractal  dimension  for  the  chosen  orientation  is  the  average  value 
of  those  obtained  fractal  dimensions  for  different  sets. 

The  average  fractal  dimension  measured  in  different  orientations  13  is  shown  in 
Table  5.6.  The  statistical  analysis  of  the  data  will  be  explained  in  Appendix  B.  From 
the  data,  it  is  seen  that  the  fractal  dimension  varies  not  only  with  the  fracUire  planes  but 
also  with  the  tensile  surfaces  under  flexure.  Thus  the  change  in  orientations  includes 
both  the  choice  of  the  fracture  plane  and  that  of  the  tensile  surface.  This  implies  that  the 
fractal  process  varies  with  the  fracture  plane  and  tensile  surface  plane,  which  further 
implies  that  the  fracmre  process  in  single  crystal  Si  is  self-affme.  Feamres  in  a  self- 
affine  object  are  related  by  a  vector  (scaling)  quantity  in  different  directions. 

The  {111}  fracture  plane  is  relatively  smooth,  while  the  {100}  is  torUious.  The 
{110}  is  somewhere  in  between.  Fig.  5.16  shows  feaUires  throughout  the  {110}  fracture 
surface.   If  a  small  area  in  the  lOOX  pictore  is  magnified  to  1,000X,  similar  feamres  are 
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(a) 


(b) 


Fig.  5.16  Self-similar  features  shown  on  the  fracture  surface  of  the  {110}  fracture 

plane  with  the  {100}  tensile  surface.  If  a  small  area  in  the  lOOX  picture  is  magnified  to 
1,000X,  the  same  features  are  seen.  When  a  small  region  in  the  1,000X  is  magnified  to 
10,000X,  once  again  the  same  features  appear.  Scale  invariance  is  displayed,  (a)  A  lOOX 
picture,  (b)  A  1 ,000X  picture,  (c)  A  10,000X  picture. 
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(c) 


Fig.  5.16  Continued. 
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seen.  When  a  small  region  in  the  1,000X  is  magnified  to  10,000X,  once  again  similar 
features  appears.  This  is  a  suggestion  of  self-similarity  and  scale  invariance.  Self- 
similarity  means  that  feamres  in  different  regions  appear  to  be  similar  to  one  another. 
They  can  be  related  by  a  scalar  multiple.  Scale  invariance  occurs  if  a  feature  on 
different  scales  appears  to  be  the  same.  These  properties  indicate  that  a  fractal  can  be 
created  by  using  a  shape  to  be  repeated  and  reduced  in  size  followed  by  a  prescribed 
sequence.  Thus,  those  photographs  of  the  fracture  surfaces  of  single  crystal  silicon 
display  fractal  patterns. 

Table  5.6    Fractal  dimension  for  different  orientation 


fracture  plane 

tensile  surface 

fractal  dimension 

{100} 

{110} 

1.16  ±0.04 

{110} 

{100} 
{110} 

1.10  ±0.04 
1.04  ±0.03 

{111} 

{110} 
{112} 

1.06  ±0.02 

1.05  +  0.00 

'                 1 

Fracture  and  fractal  process  varies  with  the  fracture  plane  and  tensile  surface 
plane.  The  fracture  process  in  single  crystal  is  self-affine  when  the  fracture  plane  or  the 
tensile  surface  are  different.  Self-similarity  and  scale  invariance  just  occurs  in  the  same 
fracture  plane  with  the  same  tensile  surface. 

If  the  fracture  surfaces  for  different  orientations  are  carefully  compared,  the 
topography  looks  different.  It  suggests  that  self-similar  behavior  just  occurs  in  the  same 
orientation.  It  is  relatively  easy  to  explain  this  phenomenon  from  the  easiest  cleavage 
plane  point  of  view.  The  fracture  path  tends  to  follow  the  least  energy  path  and  thus  the 
crack  will  try  to  find  the  easiest  cleavage  plane  at  the  time  of  propagation  or  choose  a 
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plane  which  consumes  the  lowest  energy.  Here,  the  easiest  cleavage  plane  might  not  be 
necessarily  the  {111}  plane  because  the  {111}  plane  is  not  always  available  during  the 
fracture  process.  The  crack  path  may  follow  a  plane  which  contains  a  resolved  stress 
high  enough  to  fracture  and  at  the  same  time  consumes  the  least  energy  after  the 
external  stress  is  applied.  Thus,  if  the  orientation  changes,  the  path  of  fracture  can 
change.  Since  the  fracture  path  or  angle  relative  to  the  designed  fracmre  plane  changes, 
the  topography  of  the  fracmre  surface  will  also  change  and  the  fractal  dimension  will  be 
different  for  those  different  orientations.  Self-similarity  is  exhibited  on  the  fracmre 
surface  for  the  same  fracture  plane  and  loading  directions,  while  self-affinity  is  exhibited 
as  the  loading  direction  changes. 

Relationship  Between  Toughness  and  Fractal  Dimension 

The  relationship  between  the  toughness  and  fractal  dimension  is  one  of  the  main 
objectives  in  this  stady.  From  Table  5.7,  we  can  see  that  the  moduli  of  the  fracmre 
planes  do  not  relate  to  the  toughnesses  as  the  orientation  changes.  Generally,  materials 
with  a  higher  modulus  are  stiff er  and  have  higher  toughness.  But  this  is  not  the  case 
here  for  single  crystal  Si.  Indeed,  the  easiest  cleavage  plane  with  the  lowest  toughness 
which  is  the  {111}  plane  possesses  the  highest  modulus.  The  same  condition  occurs 
with  Ge.  Thus,  modulus  is  not  the  only  determinant  factor  for  the  toughness  in  single 
crystal  silicon  or  anisotropic  materials. 

From  Table  5.7,  we  can  see  that  the  fractal  dimension  increases  as  the  critical 
stress  intensity  factor  increases,  which  means  that  a  fracmre  surface  possessing  higher 
tormosity  indicates  that  the  fracmre  produces  more  surface  during  crack  propagation. 
Here,  an  example  is  given  to  illustrate  this  point  of  view.  If  there  is  a  triangle 
originally  with  area,  A,  and  perimeter  P  and  each  side  is  replaced  by  one  specific 
generator  as  shown  in  Fig.  5.17,  the  perimeter  and  area  of  the  new  object  would  be 
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(a) 


is  replaced  by 


Area  =  A 
Perimeter  =  P 


(b)        After  the  first  generation 


After  the  second  generation 


Area  =  A(  1  +  1/3) 
Perimeter  -  P  (4/3) 


Area  =  A  (1  +  1/3 +  (1/3)2) 
Perimeter  =  P  (4/3)2 


Fig.  5 . 1 7  Construction  of  a  fractal  obj  ect. 
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(c)  After  nth  construction 


Pnth     =Px(4/3)n 

=  infinity,  if  n  is  a  very  large  number; 

Anth    =  A  X  ( 1  +  1/3  +  (1/3)2  +  (i/3)3  +  ...  +  (i/3)n ) 

=  A  X  ( 1 .5  ),         if  n  is  a  very  large  number. 


Fig.  5.17 


Continued. 
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equal  to  P  x  (4/3)  and  A  x  (1  +  1/3)  separately  after  the  first  generation.  The 
perimeter  and  area  of  the  new  object  would  be  equal  to  P  x  (4/3)  and  A  x  {1  +  1/3  + 
(1/3)^},  separately  after  the  second  generation.  After  the  nth  generation,  the  new 
perimeter  would  become 

Pnth     =Px(4/3)n-l  (5.8) 

->  infinity,  if  n  is  a  very  large  number; 

while  the  area  of  the  generated  object  would  become 

Anth     =  A  X  (  1  +  1/3  +  (1/3)2  +  (1/3)3  +  ...  +  (l/3)n-l  )  (5.9) 

»  A  X  ( 1.5  ),  if  n  is  a  very  large  number. 

Thus,  the  more  tortuous  the  object,  the  more  surface  area  generated.  Using  this 
construction,  an  object  with  finite  area  but  infinite  perimeter  could  be  generated  and  the 
fractal  dimension  of  this  object  is  1.26  (  =  log  4  /  log  3  as  described  in  the  fractal 
dimension  section  in  Chapter  2). 

Table  5.7    Modulus,  toughness,  surface  energy,  and  fractal  dimension  in 
different  orientations  for  single  crystal  silicon. 


fracture 
plane  (hkl) 

tensile 
surface 

Ehkl 
(GPa) 

K 

(MPa  m^^^) 

Y 
(J/m^) 

(average) 

{110} 

{110} 

169 

1.07 

3.39 

0.04 

{111} 

{112} 

188 

1.16 

3.58 

0.05 

{111} 

{110} 

188 

1.17 

3.64 

0.06 

{110} 

{100} 

169 

1.23 

4.48 

0.10 

{100} 

{110} 

130 

1.26 

6.11 

0.16 

If  this  construction  of  an  object  between  one  and  two  dimensions  is  extended  to 
two  and  three  dimensions,  a  solid  body  with  infinite  surface  could  be  constructed. 
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However,  if  we  limit  the  scale  of  generation  of  fractal  objects  to  atomic  size,  i.e., 
about  several  Angstroms,  the  steps  of  the  generation,  n,  would  not  be  very  large  and 
the  fracture  surface  with  higher  fractal  dimension  would  have  higher  surface  area  than 
one  with  lower  fractal  dimension.  Here,  this  idea  is  applied  to  the  generation  of 
fracture  surfaces.  Some  fracmre  surfaces  are  smoother  than  others,  like  the  {111} 
fracture  plane  is  smoother  (D  =  1.05)  than  the  {110}  (D  =  1.08,  the  higher  of  the  two 
chosen  tensile  surfaces)  and  {100}  (D  =  1.16,  the  higher  of  the  two  chosen  tensile 
surfaces)  fracture  plane  in  single  crystal  silicon.  In  single  crystal  Si,  there  is  one  kind 
of  bond  and  the  energy  of  the  covalent  bond  is  the  same  for  different  orientations. 
More  surface  area  means  more  bonds  are  broken  and  more  energy  is  consumed  during 
the  fracture  process.  Thus,  it  is  reasonable  to  say  that  the  {111}  plane  is  the  fracture 
plane  with  the  least  toughness  since  it  has  the  smoothest  fracture  surface  with  fewest 
bonds  broken,  while  {110}  and  {100}  planes  are  tougher  since  then:  fracture  surfaces 
show  more  irregularity  and  have  more  bonds  broken. 

Note  that  the  area  beyond  crack  branching  should  not  be  considered  in  the 
fractal  analysis  at  the  scale  for  which  we  examine  the  surface.  Thus,  only  the  fracture 
surfaces  generated  prior  to  crack  branching  were  considered  and  compared  with  the 
tortuosity  or  fractal  dimensions.  If  we  said  that  the  fracture  surfaces  on  the  {100} 
plane  are  more  tortoous  than  those  on  the  {111}  plane,  then  the  fracture  surfaces  to 
which  we  refer  are  in  the  area  before  crack  branching.  The  areas  beyond  crack 
branching  always  look  very  smooth,  because  the  "mirror"  region  repeats  itself  after 
branching,  because  of  the  self-similar  nature  of  the  fracmre  process. 

In  order  to  prevent  getting  an  infinite  value  of  ao  calculated  from  Eq.  (2.22)  if 
D*  approaches  zero,  we  modify  Eq.  22: 

y  =  Yo  +  aoEhkiD*/2,  (5.10) 

where  y  is  the  fracture  energy, 

Yo  is  the  surface  energy  for  an  Euclidean  fracture  surface  (D*=0), 
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Eq  is  a  characteristic  (atomic)  length,  and 

EjiIjI  is  the  Young's  modulus  of  the  hkl  fracture  plane. 

The  plot  of  y  vs.  ED*/2  for  different  orientations  of  single  crystal  silicon  is 
shown  in  Fig.  5.18.  The  slope  with  a  constant  value  suggests  that  a.^  is  a  constant, 
which  means  that  fracture  surface  energy  is  the  sum  of  the  energy  necessary  to  break 
bonds  and  the  energy  which  makes  the  flat  surface  into  a  tortuous  surface.  The  latter 
energy  is  consumed  in  the  fracture  process  and  contributes  to  the  formation  of  a 
fracture  surface  geometry.  This  meaning  agrees  very  well  with  the  fact  that  the 
fracture  surface  is  never  ideally  smooth  even  in  the  easiest  cleavage  plane.  By  using 
Eq.  (5.10),  the  a^  value  is  found  to  be  4.1  A,  which  is  interpreted  to  be  the  stretched 
length  of  the  Si-Si  bond  just  before  fracture. 

As  seen  in  Fig.  5.18,  the  value  of  y^  in  Eq.  (5.10)  is  the  intersection  of  the  line 
with  the  Y  coordinate  at  ED*  =  0  and  is  found  to  be  1.7  J/m^,  which  is  very  close  to 
the  value  of  surface  energy  calculated  from  bond  breaking  as  discussed  in  the  fracture 
toughness  measurement  section.  Those  calculated  values  are  found  to  be  2.3  J/nfi  for 
the  {100}  plane,  1.67  J/m2  for  the  {110}  plane,  and  1.36  J/m2  for  the  {111}  plane. 
The  assumption  used  to  calculate  the  surface  energy  in  the  fracture  toughness  section  is 
that  the  bonds  all  break  in  one  plane  with  no  crack  branching  and  no  non-coplanar 
crack  growth.  Thus,  the  calculated  value  should  be  lower  than  the  measured  value 
from  experimental  work  and  should  be  close  to  the  energy,  y^,  needed  to  produce  a 
fracture  surface  with  Euclidean  geometry. 
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Fig.  5.18  The  plot  of  y  vs.  ED*/2.   The  slope  of  the  fitted  line  is  found  to  be  4.lA. 

The  line  intersects  with  the  y  coordinate  at  Jq,  1 .70  J/m  . 
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Results  From  Molecular  Dynamics  Simulation 

Validity  of  the  Applied  Potential 

In  this  study,  before  we  assume  that  molecular  dynamics  (MD)  is  an  applicable 
method  to  study  atomic  fracture  behavior  of  single  crystal  silicon,  the  accuracy  and 
appropriateness  of  the  numerical  procedure  and  the  interpretation  of  the  simulation 
results  should  be  discussed.  Here  we  concentrate  on  the  study  of  the  fracture  process 
at  the  atomic  level  and  the  fractal  analysis  of  the  generated  fracture  surfaces.  Earlier 
works^"^^'^^  in  MD  simulation  of  fracture  concentrated  on  silicate  glasses.  Soules  and 
Busby^^  used  molecular  dynamics  to  study  the  rheological  properties  and  fracmre 
behavior  of  a  sodium  silicate  glass.  The  work  was  followed  by  works  of  Soules, 
Ochoa  and  Simmons,^  and  Ochoa  et  al,^^  which  demonstrated  the  importance  of 
dynamic  effects  in  the  fracture  of  silicates.  A  work  by  Simmons  et  al.^^  examined  the 
mechanisms  involved  in  an  individual  bond  fracture  to  describe  a  process  by  which  an 
entke  fracture  surface  may  be  stabilized.  The  MD  technique  has  been  used  to  perform 
the  study  of  fracture  behavior  and  should  be  useful  for  the  study  of  silicon.  However, 
the  previous  studies  only  studied  ionically  bonded  materials.  This  study  includes 
covalently  bonded  materials.  Thus,  a  different  potential  had  to  be  developed  and 
Stillinger  and  Weber24  suggested  a  potential  to  be  used  for  silicon. 

Kieffer  and  Angell^^  proposed  a  way  to  relate  the  fractal  dimension  with  the 
density  function  of  aerogel  materiel  using  the  MD  technique.  They  claimed  that  a  non- 
integer  dimension  is  a  characteristic  feature  of  the  aerogel  strucmre  used  in  their  study. 
The  fractal  dimension  can  be  associated  with  dynamic  properties  in  relation  to  many 
kinetic  processes  and  rate  phenomena.  However,  because  of  the  way  in  which  they 
obtained  the  fractal  dimension  by  using  the  change  of  the  slope  in  the  radial  distribution 
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functions  with  respect  to  the  density,  their  results  were  not  reasonable.  Their  results 
show  that  a  denser  aerogel  has  a  lower  fractal  dimension.  This  results  contradicts  the 
fact  in  nature  that  a  fully  occupied  volume  has  a  dimension  3  while  a  volume  with 
porosity  has  a  dimension  less  than  3.  Since  it  was  concluded  that  their  approach  was 
incorrect,  we  decided  to  develop  a  different  approach  from  theirs  to  obtained  the  fractal 
dimension  in  this  study. 

As  discussed  in  Chapter  3,  a  modification  of  the  Stillinger-Weber  (SW) 
potential  is  needed  to  correctly  predict  fracture  behavior  in  covalently  bonded 
materials.  Before  fracture,  silicon  possesses  a  DC  structure  which  can  be  simulated 
using  the  SW  potential  which  performs  very  well  in  covalently  bonded  and  tetrahedral 
structures.  However,  after  fracture,  the  local  structure  changes  and  those  Si  atoms 
around  the  crack  do  not  have  4  neighbors  anymore  and  some  electric  charge  must  be 
added  to  those  atoms.  A  different  potential  must  be  introduced  after  fracture  in  order 
to  describe  the  local  changes.  A  modified  Born-Mayer  potential  works  very  well  in 
silica  crystal  and  is  used  here.  Thus,  the  SW  potential  and  the  modified  Born-Mayer 
potential  were  applied  in  the  system  for  this  study.  In  Chapter  3  the  potential  and  the 
detailed  simulation  procedures  have  been  described.  The  two-body  potential  and  force 
curves  for  single  crystal  silicon  are  shown  in  Fig.  5.19.  From  the  curves,  it  is  found 
that  the  equilibrium  position  for  Si  is  about  2.34  A,  and  the  distance  for  the  cut-off 
potential  is  3.78  A.  Fig.  5.20  shows  the  magnimde  of  the  Coulombic  modified  Born- 
Mayer  potential  related  to  the  two-body  part  of  the  SW  potential.  Notice  that  the 
magnitudes  of  the  two  potentials  used  in  this  study  are  of  the  same  order. 

The  validity  of  the  numerical  procedure  and  the  interpretation  of  the  simulation 
results  can  be  checked  by  obtaining  physical  properties  from  the  MD  results  and 
comparing  these  with  the  experimental  results.  The  SW  potential  has  been  accepted 
widely  in  several  applications  such  as  thermodynamics  properties, 24,56  phase 
diagrams, 56  and  defect  study. ^^  xhe  validity  of  using  this  potential  (the  combination  of 
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Fig.  5.19  Potential  energy  and  force  curves  for  the  Stillinger- Weber  potential. 
From  the  curve,  it  is  found  that  the  equilibrium  position  for  Si  is  about  0.34 A,  and  the 
distance  for  cut-off  potential  is  3. 78 A. 
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Force 


Stillinger- Weber  force 

Coulombic  force  1-1 

Coulombic  force  2-1 

Coulombic  force  2-2 

Coulombic  force  3-1 

Coulombic  force  3-2 

Coulombic  force  3-3 

Coulombic  force  4-1 

Coulombic  force  4-2 


Fig.  5.20  The  magnitude  of  the  Coulombic  modified  Born-Mayer  potential  is 
related  to  the  two-body  part  of  the  SW  potential.  It  is  found  the  magnitudes  of  the  two 
potentials  used  in  this  study  are  of  the  same  order. 
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the  SW  and  the  Coulombic  modified  Born-Mayer  potential)  as  applied  to  single  crystal 
silicon  in  this  study  will  be  checked  later  with  the  pair  correlation  functions  and  system 
pressure  at  p  =  2.33  gm/cm^.  Generally,  before  the  fracture  process  takes  place,  only 
the  SW  potential  interacts  between  atoms  because  each  Si  atom  has  four  neighbors. 
Thus  in  doing  evaluations  of  the  system  pressure  and  pair  correlation  functions,  the  SW 
potential  is  the  main  factor  at  low  temperature  and  before  fracture. 

Since  the  judgment  of  appropriate  potentials  rests  on  the  fit  with  experiment,  we 
compared  the  pair  correlation  function  obtained  using  the  SW  and  modified  Born- 
Mayer  potentials  with  measured  pair  correlation  functions.  Fig.  5.21  shows  the  pair 
correlation  function  of  single  crystal  Si  at  room  temperature  from  MD  simulation.  Fig. 
5.22  shows  the  same  (simulated)  function  at  a  temperature  a  little  lower  than  the 
melting  point.  Fig.  5.23  shows  the  same  function  at  a  temperature  a  little  over  the 
melting  point.  In  Fig.  5.22,  the  curve  from  the  simulation  fits  quite  well  with  the 
measured  one  at  the  first  and  second  peaks. ■^'*  After  the  second  peak,  some  differences 
show  up  because  only  two-body  and  three-body  potentials  were  considered  in  the  SW 
potential  and  thus,  it  is  not  really  a  long-range  potential.  So  the  agreement  at  the  first 
and  second  peaks  is  sufficient  for  local  structure  calculations. 

Fig.  5.24  shows  the  pressure  of  the  system  at  room  temperature  if  the  actual 
density  of  silicon  is  selected  in  the  input-data  file.  The  system  remains  balanced  slightly 
above  0  GPa  with  some  flucmations  caused  by  thermal  vibrations.  Cohesive  energy, 
crystal  to  liquid  phase  transformation,  and  structure  factors  obtained  from  MD 
simulations  have  been  studied  to  compare  with  the  experimental  values  by  Stillinger 
and  Weber^^  and  are  found  to  fall  in  a  reasonable  range. 
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Fig.  5.21         Pair  correlation  function  of  silicon  at  room  temperature  from  MD 
simulation. 
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Fig.  5.22         Pair  correlation  function  for  single  crystal  silicon  at  a  temperature 
slightly  lower  than  the  melting  point. 
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Fig.  5.23         Pair  correlation  function  for  single  crystal  silicon  at  a  temperature 
slightly  higher  over  the  melting  point. 
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Fig.  5.24         Pressure  of  the  unstrained  system  at  room  temperature  at  p  =  2.33. 
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Fracture  Using  MP  Simulation 

Analysis  of  fracture  by  MD  is  performed  by  recording  system  properties  and 
structural-related  data  during  expansion.  The  primary  interest  in  the  system  properties 
is  focused  on  the  system  pressure.  The  pressure  data  can  be  used  to  determine  the 
stress-strain  curve.  During  simulations,  the  location  of  fracture  varied  for  different 
simulations  even  when  the  conditions,  like  strain  rate,  number  of  atoms,  temperature, 
thermal  properties,  and  density,  are  set  to  be  the  same.  However,  the  stress-strain 
curve  obtained  is  almost  exactly  the  same  for  different  simulations  with  the  same 
conditions.  This  means  that  although  the  locations  for  fracture  may  vary  due  to  random 
thermal  vibrations,  statistically  the  whole  system's  properties  will  be  the  same  despite 
the  local  differences. 

In  order  to  demonstrate  the  use  of  these  tools  in  the  investigation  of  fracture  by 
MD,  an  example  of  a  system  at  0.2  ps  strain  rate  along  the  x-axis  was  examined. 
The  system  pressure  during  0.2  ps'^  strain  rate  pulling  along  the  x-axis  is  shown  in 
Fig.  5.25.  The  system  pressure  increases  initially  then  decreases  after  reaching  some 
point.  This  indicates  that  the  tensile  stress  is  developed  in  the  system  as  the  system  is 
elongated,  as  would  be  expected.  Note  that  the  system  pressure  drops  when  the 
elongation  reaches  45  % .  This  indicates  that  the  material  has  some  kind  of  structural 
transition.  The  maximum  tensile  pressure  attained  by  the  material  is  assumed  to  be  the 
strength  of  the  material  and  at  that  moment  the  material  has  begun  to  fracture. 
However,  the  material  does  not  fracmre  in  one  instant  initially,  the  bonds  break 
primarily  individually  and  the  pressure  decreases  as  the  number  of  broken  bonds 
increases.  The  next  generation  of  broken  bonds  will  occur  near  the  already  broken 
bonds  to  form  the  crack  surface. 
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Fig.  5.25  The  system  pressure  (stress)  during  0.2  ps      strain  rate  pulling  along 

the  X-axis. 
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The  structural  changes  during  elongation  and  fracture  can  also  be  studied  from 
the  simulated  data.  The  first  peak  in  the  pair  correlation  functions  all  shift  to  greater 
distances  during  expansion  of  the  sample,  as  shown  in  Fig.  5.26. 

The  Si-Si-Si  bond  angle  distributions  can  be  obtained  from  simulations  and  are 
shown  in  Fig.  5.27.  Initially  the  distribution  has  a  peak  at  about  109°  in  agreement 
with  the  ideal  DC  structore.  After  uniformly  pulling,  the  distribution  undergoes  a 
change  which  depends  on  the  strain  rate. 

The  structure-related  data  recorded  during  simulations  are  like  frame-by-frame 
pictures  of  the  fracture  process  at  the  atomic  level.  The  MD  simulations  used  in  this 
work  are  based  on  a  program  written  by  Swiler.^^  From  these  data,  pair  correlation 
functions,  bond  strain  distributions  and  a  3-D  presentation  of  the  structure  can  be 
obtained.  Stereo  pairs  of  video  frozen  frames  give  a  visual  presentation  of  the  fracmre 
process.  The  location  of  the  atomic  sites  in  an  unstrained  single  crystal  silicon  in  the 
{100}  fracture  plane  is  shown  in  Fig.  5.28.  At  45%  strain,  fracture  is  in  progress  as 
evidenced  by  the  presence  of  a  free  volume  sphere  in  the  center  region  as  shown  in 
Fig.  5.29.  At  50%  elongation,  the  fractured  region  increases  as  evidenced  by  the 
presence  of  the  large  free  volume  and  several  broken  bonds  as  shown  in  Fig.  5.30. 
Thus,  the  ability  of  the  material  to  tolerate  stress  is  not  immediately  lost  when  fracture 
occurs  in  part  of  the  specimen  during  the  MD  simulation.  During  the  fracture  event  in 
MD  simulation,  it  has  been  observed  that  atom  particles  are  emitted  during  fracture. 
This  phenomenon  agrees  with  the  experimental  work  done  by  K'Singam  et  al.'^^ 

I  have  just  presented  a  brief  discussion  of  MD  analysis.  I  will  next  present  the 
values  of  toughness  and  fractal  dimensions  obtained  from  the  simulations  for  different 
orientations.  The  strength  dependence  on  the  strain  rate,  and  the  strength  dependence 
on  the  crack  size  will  be  discussed  later. 
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Fig.  5.26         Pair  correlation  function  during  extension.     The  first  peak  in  the  pair 
correlation  functions  all  shift  to  greater  distances  during  expansion  of  the  sample. 
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Fig.  5.27  The  Si-Si-Si  bond  angle  distribution  during  pulling.  Initially  the 
distribution  has  a  peak  at  about  109°.  After  uniformly  pulling,  the  distribution 
undergoes  a  change  which  depends  on  the  applied  strain. 
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Fig.  5.28         Unstrained  single  crystal  silicon  which  is  designed  to  fracture  in  the 
{100}  plane. 


126 


Fig.  5.29         At  47%  extension,  fracture  occurs.     The  evidence  for  fracture  is  the 
presence  of  a  free  volume  sphere  in  the  center  region. 
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Fig.  5.30         At  50%   elongation,  fracture  continues.      The  free  volume  spherical 
region  gets  larger  and  bonds  break. 
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Toughness  in  Different  Orientations 

Stress-strain  curves  as  shown  in  Fig.  5.31  have  been  obtained  in  three  different 
orientations,  i.e.,  {100},  {110}  and  {111}.  From  the  slope  of  the  stress-strain  curve, 
the  modulus  is  found  to  be  92  GPa  for  {100},  130  GPa  for  {110}  and  138  GPa  for 
{HI}.  For  single  crystal  Si,  s^,  Si2,  and  S44  are  measured  to  be  0.768,  -0.214,  and 
1.256  X  10-11  Pa-17^  Thus,  the  experimental  E^ki  values  are  calculated  to  be  130  GPa 
for  {100},  169  GPa  for  {110}  and  188  GPa  for  {111},  respectively  (Appendix  C).  The 
simulated  values  are  about  three  quarters  of  the  measured  values.  Although  the 
simulated  values  are  not  the  same  as  the  actual  values,  the  ratio  of  the  values  are  about 
the  same. 

The  fracture  strength  from  MD  simulation  is  38  GPa  for  {100},  36  GPa  for 
{110}  and  30  GPa  for  {111}  if  a  fixed  strain  rate,  0.2  1/ps,  is  applied  for  each 
orientation.  Generally,  the  strength  for  ideal  materials  is  about  E/27i.  The  simulated 
value  is  approximately  between  E/4  to  E/6  if  the  actual  modulus  is  used  or  is  about 
between  E/3  to  E/4. 5  if  the  simulated  modulus  is  used.  The  reason  for  such  high 
fracmre  strengths  in  the  sunulation  is  that  the  Stillinger-Weber  potential  has  a  3.76  A 
cutoff  distance.  Since  the  cutoff  distance  is  so  large,  the  strain  at  the  time  of  Si 
fracture  appears  to  be  unreasonably  high  for  brittle  materials  such  as  Si.  However, 
local  stresses  are  expected  to  be  near  the  theoretical  values,  which  are  high.  The 
simulation  is  considered  to  model  local  stresses. 

The  toughness  can  be  estimated  from  the  stress-strain  curve,  if  we  assume  that 
the  area  under  the  stress-strain  curve  at  the  point  where  it  reaches  the  maximum  stress 
is  the  energy  consumed  to  fracmre  the  samples.  Thus,  the  equation  to  calculate  the 
toughness  is  as  follows: 

toughness  =   J^max^j  ds.  (5.10) 
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Fig.  5.31         Stress-strain  curve  for  three  different  orientations,   {100},   {110}  and 
{111}. 
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It  is  easily  seen  from  Fig.  5.31  that  the  {100}  plane  has  the  highest  value,  while 
the  {111}  plane  has  the  lowest  value.  The  values  of  toughness  are  found  to  be  0.90d/, 
0.6 Id/,  and  0.43d/  for  the  {100},  {110}  and  {111}  planes,  respectively,  where  d/  is  the 
unit  length.  The  choice  of  the  unit  length  determines  the  value  of  the  actual  toughness 
and  will  be  illustrated  after  the  fractal  analysis  has  been  discussed. 

Since  the  simulated  moduli  are  lower  than  the  actual  values,  it  is  reasonable  to 
assume  that  the  fracmre  strength  and  toughness  will  be  lower  than  the  actual  values 
also.  However,  the  simulated  strength  is  a  little  higher  than  the  actual  value  and  the 
toughness  is  about  the  same,  if  a  reasonable  unit  length,  d/,  is  chosen.  Thus,  the  actual 
values  for  the  simulated  moduli  do  not  appear  to  be  as  important  as  other  variables 
such  as  the  (SW)  potential  truncation  and  the  two  potential  approach. 

Strength  Dependence  on  the  Strain  Rate 

The  dependence  of  the  strength  on  the  strain-rate  was  investigated  for  the  {100} 
fracture  plane.  The  stress-strain  curves  for  strain  rates  of  0.1,  0.2,  0.5,  1,2,  and  5  ps" 
^  are  shown  in  Fig.  5.32.  As  the  strain  rate  increases,  the  strength  increases  in 
agreement  with  the  experimental  work  of  Ebrahimi.76  There  are  two  reasons  for  this. 
When  specimens  are  elongated  at  the  higher  strain  rate,  more  cracks  initiate  and  several 
fracture  planes  are  formed  after  final  fracture.  So,  at  higher  strain  rates,  more  energy 
has  to  be  consumed  in  order  to  generate  more  cracks  and  fracture  surfaces.  Another 
reason  is  that  structural  relaxation  is  allowed  to  occur  at  lower  strain  rates,  but  not  at 
higher  strain  rates.  This  relaxation  allows  for  the  relief  of  stress.  Thus,  the  selection 
of  strain  rates  relative  to  the  time  scale  of  the  material  is  critical  to  the  type  of  behavior 
that  is  observed. 

A  useful  unit  of  time  at  the  atomic  scale  is  the  period  of  thermal  vibration, 
which  is  an  atom's  response  time  to  a  change  in  the  environment.     The  thermal 
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Fig.  5.32         The  strength  dependence  on  the  strain-rate.    Strength  increases  as  the 
strain  rate  increases. 
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vibration  period  is  found  to  be  0.0766  ps  for  Si.  Thus  a  0.2  ps"^  strain  rate  translates 
to  a  0.0152  expansion  rate  (unitless)  in  terms  of  the  vibration  periods,  and  a  2  ps" 
strain  rate  translates  to  a  0.152  expansion  rate.  Those  numbers  mean  that  for  the 
slower  strain  rate,  the  sample  is  elongated  1.52%  per  vibration  period,  while  at  the 
higher  strain  rate,  the  sample  is  elongated  15.2%  per  vibration  period.  Although  these 
strain  rates  sound  fast,  relaxation  by  thermal  vibrations  are  allowed  at  the  lower  strain 
rates  because  the  elongation  per  vibration  period  is  relatively  quite  small,  while  for  the 
higher  strain  rate,  the  time  for  relaxation  is  not  long  enough.  If  each  tensile 
deformation  is  small  enough  for  the  system  to  relax  before  the  next  extension,  the 
system  can  reach  its  equilibrium  condition  during  the  period  between  two  "loadings." 
Reaching  the  equilibrium  condition  means  that  the  lowest  energy  state  and,  hence,  the 
stress  state  for  the  system  at  that  moment  is  relieved.  Thus,  a  system  with  a  slower 
strain  rate  will  produce  a  lower  fracture  strength. 

After  reaching  the  maximum  stress,  the  stress  in  the  stress-strain  curve  drops 
faster  for  the  slow  strain  rate  than  the  fast  strain  rate.  The  reason  for  this  rapid  drop  is 
that  the  distance  increment  at  each  time  step  for  the  slower  strain  rate  is  small  enough 
for  the  system  to  reach  equilibrium  at  a  shorter  time.  In  the  fast  strain  rate  case,  the 
distance  increment  between  two  atoms  at  each  extension  increment  is  so  large  that  the 
system  caimot  respond  fast  enough.  Thus,  the  stress  drops  faster  if  the  strain  rate  is 
slower.  "76 

The  structural  changes  during  elongation  and  fracture  can  be  studied  for  the 
above  two  cases,  i.e.,  for  the  0.2  ps"^  and  the  2  ps'^  strain  rate.  The  first  peaks  for 
10%,  20%,  30%  and  40%  strain  in  the  pair  correlation  functions  for  the  0.2  ps'^  strain 
rate  all  shift  to  greater  distances  and  there  is  no  big  shape  change  during  expansions  of 
the  sample,  as  shown  in  Fig.  5.33.  However,  the  shape  and  position  of  the  second  and 
third  peaks  change  greatly  as  compared  to  the  first  peak.  This  large  change  means  that 
the  long  range  order  has  changed  greatly.    As  for  the  condition  with  a  2  ps'^  strain 
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Fig.  5.33         The  pair-correlation  function  for  a  0.2  ps'^strain  rate.    The  first  peak  in 
the  pair  correlation  functions  shift  to  greater  distances  during  expansion  of  the  samples. 
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rate,  as  shown  in  Fig.  5.34,  the  position  of  the  first  peak  only  changes  very  little  but 
the  width  of  the  first  peak  changes  very  much.  In  addition,  the  tip  of  the  first  peak  is 
very  sharp,  which  means  that  the  expansion  rate  is  too  fast  for  the  system  to  smooth  the 
first  peak.    The  shape  and  position  of  the  second  and  third  peaks  change  very  much 

also. 

The  Si-Si-Si  bond  angle  distribution  for  the  0.2  ps"^  and  the  2  ps'^  strain  rates 
can  be  obtained  and  are  shown  in  Figs.  5.35  and  36.  Initially  the  distribution  has  a 
peak  at  about  109°.  After  uniform  expansion,  the  shape  of  the  bond  angle  distribution 
undergoes  a  change  which  depends  on  the  applied  strain.  From  Figs.  5.35  and  36,  we 
can  see  that  the  peak  of  the  bond  angle  distribution  is  at  109°  at  first  then  it  resolves 
into  two  peaks  at  higher  strain.  While  comparing  the  shape  of  the  peak,  those  from  the 
0.2  ps"^  are  much  smoother  than  those  from  2  ps"^  The  system  with  the  slower  strain 
rate  can  reach  its  equilibrium  faster  to  produce  a  smoother  curve  than  that  with  higher 
strain  rates.  Since  the  system  with  the  slow  strain  rate  can  reach  its  equilibrium  in  a 
more  realistic  manner  than  the  faster  strain  rates,  the  simulations  used  for  determining 
the  strength  dependence  on  the  crack  size,  and  toughness  were  all  performed  at  0.2  ps' 
strain  rate. 

Strength  Dependence  on  the  Crack  Size 

The  Griffith  criterion  implies  that  if  the  material  possesses  larger  cracks,  it 
would  fracmre  at  a  lower  stress  level.  Here  a  controlled  crack  will  be  introduced  in  the 
specimen.  Specimens  with  1,  2,  4,  8,  and  16  atoms  removed,  respectively,  will  be 
added  to  the  simulation  criteria.  Two  types  of  cracks  are  introduced  here,  one  is  void- 
like, and  another  is  plane-like. 

The  stress-strain  curves  for  specimens  with  a  void-like  crack  are  shown  in  Fig. 
5.37.   It  is  seen  that  the  fracture  strength  varies  with  the  size  of  the  cracks.    The  crack 
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Fig.  5.34  The  pair-correlation  function  for  a  2  ps'^ strain  rate.  The  first  peak  in 
the  pair  correlation  functions  all  shift  to  greater  distances  during  expansion  of  the 
samples. 
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Fig  5  35  The  Si-Si-Si  bond  angle  distribution  during  pulling  for  the  0.2  ps"  strain 
rate.  Initially  the  distribution  has  a  peak  at  about  109°.  After  uniformly  pulling,  the 
distribution  undergoes  a  change  which  depends  on  the  applied  strain. 
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Fig.  5.36  The  Si-Si-Si  bond  angle  distribution  during  pulling  for  the  2  ps'^  strain 
rate.  Initially  the  distribution  has  a  peak  at  about  109°.  After  uniformly  pulling,  the 
distribution  undergoes  a  change  which  depends  on  the  applied  strain. 
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Strength  dependence  on  the  crack  size.    The  shape  of  the  crack  is  void- 
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sizes,  fracture  strengths  and  calculated  critical  stress  intensity  factors  are  listed  in  Table 
5.8. 

Table  5.8  Void  sizes  and  fracture  strengths  from  MD 


#  of  atoms 
removed 

flaw  size* 
(A) 

fracture  stress 
(GPa) 

apparent  stress 

intensity  factor** 

(MPa  m^'^) 

0 

2 

38 

0.7 

1 

3 

29 

0.6 

2 

4 

27 

0.7 

4 

5 

24 

0.7 

8 

6 

22 

0.7 

16 

7 

20 

0.7 

*Flaw  size  is  considered  as  the  diameter  of  the  void  flaw,  2c. 
**Stress  intensity  factor  is  calculated  using    Kc  =  1.24  •  af  ■  c^^^. 


The  stress-strain  curves  for  the  specimens  with  a  planar  sharp  crack  are  shown 
in  Fig.  5.38.  It  is  seen  that  the  fracture  strength  varies  with  the  size  of  the  cracks. 
The  crack  sizes,  fracmre  strengths  and  calculated  critical  stress  intensity  factors  are 
listed  in  Table  5.9. 

Although  a  periodic  boundary  condition  is  applied  during  the  sunulation,  the 
simulated  results  agree  quite  well  with  linear  elastic  fracmre  mechanics  (similar  to  the 
Griffith  criterion)  for  void-like  crack  or  plane-like  crack  (sharp  crack).  The  strength  of 
the  specimen  decreases  as  the  size  of  the  crack  increases  and  the  calculated  critical 
stress  intensity  factors  for  both  cases  are  approximately  constant.  Although  the  stress 
intensity  factor  is  a  value  to  describe  the  toughness  of  continuum  solids,  it  is  applied 
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Fig.  5.38         Strength  dependence  on  the  crack  size.   The  shape  of  the  crack  is  planar 
and  simulates  a  sharp  crack. 
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here  if  the  simulated  systems  are  considered  a  combination  of  a  continuum  solid  away 
from  the  crack  region  and  individual  atoms  around  the  crack  region.  The  absolute 
value  of  the  calculated  stress  intensity  factor  is  not  really  important  here.  The 
important  observation  here  is  that  a  constant  value  of  critical  stress  intensity  factor  is 
obtained  using  MD  techniques  for  different  crack  sizes  with  all  other  parameters  the 
same. 

Table  5.9    Crack  size  and  fracture  strength  for  planar  cracks  from  MD 


#  of  atoms  removed 


flaw  size* 
(A) 


fracture  stress 
(GPa) 


apparent  stress 

intensity  factor** 

(MPa  m^''^) 


38 


OJ 


29 


0.8 


26 


0.8 


23 


0.7 


11 


20 


0.7 


12 


13 


18 


0.7 


18 


14 


16 


0.7 


^law  size,  2c,  is  the  diagonal  length  of  the  removed  rectangular  or  square  plane  of 
atoms. 


**Stress  intensity  factor  is  calculated  using  K^  =— ctVtic  s  1.13aVc  .''^  which  is  used 


for  the  embedded  penny-shaped  crack  condition. 
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Fractal  Analysis  Using  MP  Results 

The  fractal  dimension  from  the  MD  simulation  can  be  obtained  as  shown  in 
Chap.  4.  After  the  fracture  has  occurred  in  the  simulation  as  shown  in  Fig.  4.8,  we 
would  try  to  obtain  the  fracture  surface  from  the  simulation.  Each  atom  is  assumed  to 
be  a  sphere  surrounded  with  an  electron  cloud  of  diameter  3.76  A  (the  potential  energy 
cut-off  distance).  The  radius  of  each  atom  is  assumed  to  be  the  radius  of  the  potential 
field.  Thus,  the  fracmre  surface  looks  like  a  surface  with  intersecting  spheres,  as 
shown  in  Fig.  4.9.  After  obtaining  the  surface,  the  slit-island  technique  as  described  in 
the  experimental  procedure  section  in  Chap.  4  is  applied  here  to  obtain  the  fractal 
dimension.  In  order  to  get  a  slit-island,  a  cut  at  constant  height  is  applied  at  the 
fracture  surface.  After  the  cut,  the  obtained  slit-island  looks  like  a  plane  with  lots  of 
circular  disks  on  it  as  shown  in  Fig.  4.10.  The  perimeter  of  the  selected  island  is 
measured  using  different  scales  as  described  before  and  the  fractal  dimension  can  be 
obtained  using  a  Richardson  plot.  The  bubble-like  fracmre  surface  and  the  slit-islands 
for  each  fracture  plane  are  shown  in  Figs.  5.39,  40,  and  41.  It  is  easy  to  see  that  in 
Fig.  5.39,  which  is  taken  from  the  {100}  fracture  plane,  the  slit-island  looks  much 
more  irregular  than  those  from  the  {110}  and  {111}  fracture  planes  as  shown  in  Figs. 
5.40  and  5.41. 

The  fractal  dimension  was  measured  with  respect  to  different  orientations,  and 
those  data  are  shown  in  Table  5.10.  From  the  data,  it  was  found  that  the  {100} 
fracture  plane  has  the  highest  fractal  dimension,  while  the  {111}  fracmre  plane  has  the 
lowest  fractal  dimension.  If  we  compare  the  fractal  dimensions  with  the  toughness 
obtained  by  calculating  the  area  under  the  stress-strain  curve  at  the  point  it  reaches  the 
maximum  stress  level,  i.e.,  the  fracture  strength,  from  MD  simulation  as  in  Table 
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Fig.  5.39  The  slit-island  obtained  from  the  simulated  fracture  on  the  {100}  fracture 
plane.  The  slit-island  looks  much  more  irregular  than  those  from  the  {110}  and  {111} 
fracture  planes  as  shown  in  Figs.  5.40  and  5.41. 
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F.  5.40 
plane. 


The  slit-island  obtained  from  the  simulated  fracture  on  the  {110}  fracture 
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Fig.  5.41         The  slit-island  obtained  from  the  simulated  fracture  on  the  {111}  fracture 
plane. 
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5.10,  then  we  see  that  the  fracture  plane  with  higher  toughness  has  a  higher  fractal 
dimension. 

Table  5.10  Fractal  Dimensions  Determined  from  MD  Simulations 


Fracture 
Plane 

(MD) 

Fractal  Dimension 

(unitless) 

(MD) 
Toughness^ 

{100} 

1.16  +  0.06 

0.90  d/ 

{110} 

1.11  ±0.05 

0.61  d/ 

{111} 

1.09  ±0.04 

0.43  d/ 



''"The  toughness  is  estimated  from  the  stress-strain  curves  from  MD  simulation,  if 
toughness  is  assumed  as  the  area  under  the  stress-strain  curve  at  the  point  where  it 

reaches  the  maximum  stress,  i.e.,  fracture  strength.   Toughness  =   J^max ^^  ^^ 


Comparison  Between  Measured  And  Simulated  Results 


Before  the  results  are  compared,  the  difference  in  the  experimental  and 
simulation  techniques  should  be  pointed  out.  In  the  experimental  procedure,  the 
specimens  were  tested  using  three-point  flexure.  Thus,  not  only  the  fracture  plane,  but 
also  the  tensile  surface  will  have  some  effect  on  the  measured  toughness  and  fractal 
dimension  value.  In  the  MD  simulation,  the  specimens  were  tested  in  tension  instead 
of  bending.  Thus,  the  stress  states  were  different  between  the  experiment  and 
simulation.  In  the  bending  test,  there  will  be  a  neutral  plane  on  the  specimen.  Thus, 
while  part  of  the  sample  is  under  tensile  stress,  the  rest  is  under  compressive  stress. 
But  in  the  tension  test,  only  tensile  stress  exists  on  the  cross-section.  Since  a  periodic 
boundary  condition  is  applied  in  the  MD  simulation,  the  only  difference  in  the  chosen 
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orientations  are  the  fracture  planes.   The  choice  of  the  orientation  of  the  tensile  surface 
does  not  affect  the  results  of  the  simulations. 

In  this  study,  the  important  factor  is  the  fracture  process  near  the  crack  tip.  In 
the  experimental  work,  the  stress  state  on  the  bending  specimen  near  the  crack  tip  is 
almost  uniform  tension.  Since  the  region  near  the  crack  tip  is  nearly  the  same  for  the 
two  conditions  of  loading,  it  is  reasonable  to  compare  the  data  between  the 
experimental  and  simulation  results  even  though  the  testing  methods  were  different. 

Another  difference  between  the  experimental  and  simulation  procedure  is  the 
size  of  the  specimen.  In  the  experimental  work,  the  size  of  the  specimen  cross-section 
varies  between  about  3mm  •  5mm  to  6mm  •  6mm.  The  scale  used  to  measure  the 
fractal  dimension  is  in  the  order  of  lam's.  But  in  the  MD  simulation,  the  size  of  the 
specmien  cross-section  is  about  20  A  •  30  A  to  50  A  ■  50  A,  which  is  about 
1/1,000,000  of  the  size  used  in  the  experimental  specimen.  The  difference  in  sizes  are 
considerable.  One  may  say  that  since  a  periodic  boundary  condition  is  applied  in  the 
MD  simulation,  then  the  size  of  the  simulated  system  is  infinite.  If  this  assumption  is 
correct  then  the  difference  in  scale  will  not  be  that  important.  The  scale  used  to 
measured  the  fractal  dimension  is  on  the  order  of  Angstroms,  which  is  about  1/1,000  of 
the  scales  used  in  the  experimental  work.  Thus,  if  the  fractal  dimensions  measured  in 
the  experimental  work  and  simulation  agree  with  each  other,  then  the  fractal  dimension 
is  the  same  at  the  jxm  level  and  at  the  A  level.  The  agreement  between  experiment  and 
simulation  implies  that  the  valid  scale  range  to  measure  fractal  dimension  is  at  least 
from  the  Angstroms  level  to  the  [im  level.  Thus,  the  fractal  phenomenon  in  the 
fracmre  process  exists  for  a  large  scale  range. 
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Comparison  of  Toughness 

The  data  obtained  from  the  experimental  and  simulation  work  is  shown  in  Table 
5.11.  The  fracture  energy,  y,  calculated  from  the  stress-strain  curve  (Fig.  5.31)  from 
the  MD  simulation  depends  on  the  d/  unit  length  chosen.  However,  the  unit  length  is 
unknown  here.  If  the  unit  length  d/  is  chosen  as  the  Si-  Si  bond  length,  2.34  A,  the 
calculated  y  is  about  half  of  the  experimental  value.  The  cutoff  distance  and  the 
characteristic  length,  a^,  are  candidates  to  be  proposed  here  as  the  unit  length  d/  in 
order  to  calculate  the  y  value.  A  comparison  of  the  calculated  values  are  listed  in  Table 

5.11.  From  those  data,  those  y  value  calculated  from  the  characteristic  length  are 
relatively  close  to  the  measured  value.  Although  we  cannot  say  the  characteristic 
length  calculated  from  the  fractal  dimension  is  definitely  the  unit  length  during  fracture, 
the  similarity  of  values  in  Table  5.11  probably  means  that  the  characteristic  length  is 
related  to  the  atomic  process  during  fracture  in  some  way. 

Comparison  Between  Fxperimental  and  Simulation  Results  for  Fractal  Dimension 

The  data  obtained  from  experimental  and  simulation  results  are  shown  in  Table 

5.12.  The  simulation  results  show  that  the  {100}  fracture  plane  has  the  highest  fractal 
dimension,  while  the  {110}  and  the  {111}  fracture  planes  have  lower  values. 

In  the  experimental  work,  the  fractal  dimensions  have  been  measured  for  two 
tensile  surfaces  in  the  {110}  and  {111}  planes.  This  made  comparison  more  difficult. 
In  the  MD  simulations,  only  tension  tests  are  performed. 
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Table  5.11       The  toughness  value  from  experimental  and  simulation  measurement 


Fracture  plane 

y  (experiment) 
(J/m2) 

y*  (MD) 
(J/m^) 

y**  (MD) 
(J/m^) 

{100} 

6 

4 

4 

{110} 

4 

3 

2 

{111} 

3 

2 

2 

*  the  unit  length,  d/,  is  the  average  characteristic  length  ao  =  4.1  A, 

calculated  from  the  measured  data. 

**the  unit  length  is  the  cutoff  distance  of  the  SW  potential  (3.8  A). 


It  is  very  clear  that  both  the  simulation  results  and  the  experimental  results  have 
the  same  trend  and  fall  in  a  reasonable  range.  The  simulated  and  experimental  values 
for  the  fractal  dimension  are  about  the  same  even  though  the  difference  in  measuring 
scale  is  large  between  the  experimental  and  simulation  results.  The  MD  technique 
provides  a  means  to  estimate  the  fractal  dimension.  Here  the  highest  dimension  was 
selected  if  two  tensile  surfaces  were  tested  for  one  fracture  plane.  Notice  that  better 
agreement  is  obtained  if  we  compare  the  column  with  the  higher  fractal  dimension  and 
the  one  obtained  from  MD  simulation  in  the  Table  5.12.  The  results  suggest  that  the 
highest  fractal  dimension  for  any  given  plane  is  the  "true"  dimension. 
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Table  5.12       Fractal  dimension  measured  from  experimental  and  simulation  work 


fracture 
plane 

tensile 
surface 

toughness 

from 

measurement 

(MpaVm) 

fractal 
dimension 

from 
experiment 

fractal 
dimension 

from 
experiment 
-  the  higher 

fractal 
dimension 

from 
simulation 

{100} 

{110} 

1.26  +  0.06 

1.16  ±0.04 

1.16 

1.16  ±0.06 

{110} 

{100} 
{110} 

1.23  ±0.08 
1.07  ±0.09 

1.10  ±0.04 
1.04  ±0.03 

1.10 

1.11  ±0.05 

{111} 

{112} 
{110} 

1.16  ±0.08 

1.17  ±0.08 

1.05  ±0.00 

1.06  ±0.02 

1.06 

1.09  ±0.04 

CHAPTER  6 
CONCLUSIONS 


This  dissertation  studied  the  fracture  of  anisotropic  single  crystal  silicon  using 
both  experimental  and  analytical  techniques.  I  first  will  summarize  the  main 
conclusions  of  the  experimental  results,  then  summarize  the  main  conclusions  of  the 
analytical  results  from  molecular  dynamics  simulations.  I  will  then  summarize  the 
comparison  of  the  experimental  and  analytical  results. 

Fractal  fracture  behavior  was  observed  in  single  crystal  silicon.  Self-similarity 
and  scale-invariance  is  suggested  by  the  fracture  surface  appearance  as  shown  in  Fig. 
5.16.  Fracture  toughnesses  and  fractal  dimensions  can  be  related  to  the  atomic  scale 
of  silicon  using  fractal  analysis.  Using  fractal  geometric  principles,  the  fractal 
dimension  on  the  fracmre  surface  allows  us  to  connect  the  macroscopic  topography  to 
the  atomic  process  from  the  calculation  of  the  characteristic  length,  a^,. 

The  fractal  dimensional  increment,  D*,  is  higher  for  the  fracture  plane  with  the 
higher  toughness  values  in  agreement  with  other  studies  which  found  higher  toughness 
values  associated  with  higher  fractal  dimensions.  ^^'^^•'*^ 

The  crack  propagation  path  for  cracks  starting  on  the  {100}  and  {110}  fracmre 
planes  both  initiate  on  the  original  plane  and  have  the  tendency  to  deviate  to  the  {111} 
plane  as  seen  in  Figs.  5.7,  5.8,  5.9  and  5.10.  This  observation  suggests  that  fracture 
would  be  more  likely  to  deviate  or  branch  on  to  those  fracture  planes  with  lower 
fracture  energy  during  crack  propagation. 

In  this  dissertation,  the  fracture  strength  dependence  on  the  crack  size,  the 
fracmre  energy  for  different  orientations,  the  fractal  dimension  obtained  from  the 
simulated  fracture  surface,  and  the  strength  dependence  on  the  strain  rate  have  been 
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modeled  using  molecular  dynamics.     For  the  first  time  a  method  to  determine  the 
fractal  dimension  on  the  MD  generated  fracmre  surface  was  proposed  here  and  found 

to  be  applicable. 

I 
The  size  of  the  introduced  crack  is  an  important  factor  for  the  simulated 

strength  in  MD  simulation.    If  a  crack  is  introduced  in  the  simulation,  the  fracture 
strength  decreases  as  the  flaw  size  increases.    The  simulated  results  agree  quite  well  f 

with  the  Griffith  criterion  for  those  specimens  with  either  void-like  flaws  or  the  plane- 
like flaws.    The  critical  stress  intensity  factor  obtained  from  simulations  with  a  plane-  j 
like  crack  (sharp  crack)  agreed  with  the  macroscopic  experimental  values. 

The  strength  increases  as  the  applied  strain  rate  increases  in  MD  simulation, 
which  agrees  with  the  experimental  work  done  by  Ebrahimi.  One  probable  reason 
for  this  behavior  is  that  more  cracks  initiate  in  the  high  strain  rate  case  and  these 
initiations  consume  more  energy.  Another  possible  reason  is  that  relaxation  is  not 
allowed  or  equilibrium  is  not  reached  in  the  high  strain  rate  case. 

The  fracture  energy  results  from  MD  simulation  agree  quite  well  with  those 
from  the  experimental  work.  The  relation  of  the  simulated  surface  energy  on  different 
planes  is  in  the  same  range  with  that  from  experimental  work.  The  present  results 
suggest  that  using  MD  simulation  to  model  fracture  behavior  is  applicable. 

The  fractal  dimension  measured  from  MD  simulation  is  about  the  same  with 
that  from  experimental  work.  It  suggests  that  using  MD  simulation  to  model  fractal 
analysis  is  acceptable.  The  MD  results  implied  that  the  fractal  dimension  can  be 
measured  using  a  scale  from  A  level  in  MD  simulation  to  fim  level  on  the  fracture 
surface. 

From  the  plot  of  y  vs.  ED*/2  for  different  orientations  in  single  crystal  silicon, 
a  constant  slope  suggests  that  ao  is  a  constant  for  single  crystal  silicon.  The 
characteristic  length,  a.^,  was  calculated  from  experimental  data  to  be  4.1  A.  This 
value  is  very  close  to  the  stretched  length  of  the  Si-Si  bond  just  before  fracture  which 
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is  estimated  from  MD  to  be  3.8A,  the  cutoff  distance  of  the  Stillinger- Weber  potential. 
From  this  work  and  previous  work,  we  know  that  fracture  energy  is  the  sum  of  the 
energy  necessary  to  breaks  bonds,  plus  the  energy  to  generate  light,  heat  and  sound, 
plus  the  energy  to  produce  a  tortuous  surface.  This  latter  energy  is  consumed  in  the 
fracmre  process  and  contributes  to  the  formation  of  a  fracmre  surface  which  is  not 
Euclidean.  The  results  shown  here  imply  that  the  fracture  surface  is  never  ideally 
smooth  even  in  the  easiest  cleavage  plane. 

The  experimental  and  simulation  data  are  summarized  in  Table  6.1. 

Table  6. 1        Summarized  results  of  this  study 


fracmre 
plane 

tensile 
surface 

(Lit.) 

E 
(GPa) 

(Exp.) 

K 

(MPa  m^'^) 

(Exp.) 

Y 
(J/m^) 

(MD) 

Y 
(J/m2) 

(Exp.) 
D* 

MD) 
D* 

{110} 

{110} 

169 

1.07 

3.39 

t 

0.04 

t 

{111} 

{112} 

188 

1.16 

3.58 

2 

0.05 

0.09 

{111} 

{110} 

188 

1.17 

3.64 

t 

0.06 

t 

{110} 

{100} 

169 

1.23 

4.48 

3 

0.10 

0.11 

{100} 

{110} 

130 

1.26 

6.11 

4 

0.16 

0.16 

'  Only  the  choice  of  the  fracmre  plane  is  important  for  those  data  from  the  MD 
simulation. 


APPENDIX  A 
SUMMARY  OF  DATA  FROM  DIFFERENT  INVESTIGATORS 


Those  fracture  toughnesses  measured  by  Oilman^,  Chen  and  Leipold,22  St. 
John,57  Jaccodine,64  Myers  &  Hillsberry,65  and  Kalwani66  are  listed  in  Table  A.l. 

Table  A.  1   Comparison  of  data  obtained  from  different  techniques  and  investigators. 


reference 

techniques 
applied 

{100} 

fracture  plane 

(MPa  ml/2) 

{110} 

fracmre  plane 

(MPa  ml/2) 

{111} 

fracmre  plane 

(MPa  ml/2) 

Jaccodine64 

double  cantilever 

0.75 

0.71 

0.62 

Chen& 
Leipold22 

indentation  & 
four  point  bend 

1.26* 

1.19* 

1.09* 

Kalwani66 

indentation 

{110}tl.25 
{100}t0.86 

{110}t0.83 
{100}t  1.01 
{111}*  0.81 

{110}*  1.22 

Oilman^ 

double  cantilever 

0.62 

Myers  & 
Hillsberry65 

notched  four 
point  bend 

1.34-2.85 

St.  John57 

Double 
cantilever 

{110}*  0.93 

this  work 

strength 
Indentation 

{110}*  1.24 

{100}*  1.19 
{110}*  1.05 

{110}*  1.19 
{112}*  1.21 

this  work 

three  point  bend 

{110}*  1.26 

{100}*  1.23 
{110}*  1.07 

{110}*  1.17 
{112}*  1.16 

*Modified  data.  (cf.  p.  90-91) 
*Indentation  plane 
*Tensile  surface. 
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APPENDIX  B 
STATISTICAL  ANALYSIS 


The  average  fractal  dimension  measured  for  different  orientations  are  shown  in 
Table.  B.l.   Those  data  are  found  to  be  very  close  to  each  other. 

Table  B.l    Fractal  dimension  for  different  orientations 


fracture  plane 

tensile  surface 

fractal  dimension 

{100} 

{110} 

1.16  +  0.04 

{110} 

{100} 
{110} 

1.10  ±0.04 
1.04  +0.03 

{111} 

{110} 
{112} 

1.06  +  0.02 
1.05  +  0.00 

For  the  {100}  fracture  plane, 


a 


the  average  fractal  dimension,  x±2ax  =  1.16 ±2^=,  and 

Vn 


y(x-x) 

the  standard  deviation,  s  =  -, — ^^ =  0.04 . 

V      n-1 


To  form  the  interval  of  2  standard  deviations  around  x ,  we  calculate 


x±2av  =  1.16±2-^«1.16  +  2    ^ 
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Here  n=  12,  Thus 


x±2av=1.16±2-^  ^1.16  +  2^  =  1.16±  0.02 


/12  V12 

Thus,  about  0.95  of  all  specunen  should  have  values  between  1.15  and  1.17. 
For  the  {110}  fracture  plane  with  a  {100}  tensile  surface, 

-    Ex 
the  average  fractal  dimension,  x  = =  1.10,  and 

n 


the  standard  deviation,  s  =  J — '^ —  =  0.04 . 

V      n-1 

To  form  the  interval  of  2  standard  deviations  around  x,  we  calculate 


x±2aj^=1.10±2-^«1.10±2    ^ 


/12  V12 

Here  n=  12,  Thus 

x  +  2ax=  1.10  +  2-^  «1.10±2^  =  1.10  +  0.02 

Vl2  Vl2 

Thus,  about  0.95  of  all  specimen  should  have  values  between  1.09  and  1.11. 

The  distribution  plots  for  these  two  cases  can  be  plotted  together  to  see  whether 
these  two  sets  of  data  are  overlapping  or  separating  as  seen  in  Fig.  B.l. 
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1.16 
s=0.04 


1.25 


Fractal  Dimension 


Fig.  B.l  The  sample  distribution  for  the  fractal  dimensions  of  the  {100}  fracture 

plane  with  a  {110}  tensile  surface,  D  =  1.16  ±  0.04,  and  the  {100}  fracture  plane 
with  a  {110}  tensile  surface,  D  =  1.10  +  0.04. 


APPENDIX  C 
YOUNG'S  MODULUS  CALCULATION 


Since  a  single  crystal  material,  in  general,  is  anisotropic,  Young's  moduli  are 
different  for  different  orientations.  The  young's  modulus  for  different  orientations  was 
calculated  from  the  stiffness  matrix  measured  by  McSkimin  and  Andreatch.'^^  single 
crystal  silicon  has  a  diamond  cubic  structure.  The  stiffness  matrix  for  a  cubic  system  is 
as  follows: 


^ij  = 


Sll  Si2  Si3  Si4  Si5  Si6 

S21  S22  S23  S24  S25  S26 

S31  S32  S33  S34  S35  S36 

S41  S42  S43  S44  S45  S46 

S5I  S52  S53  S54  S55  S56 

^61  ^62  S63  S54  855  S66 


Sll 

S12 

S12 

0 

0 

0 

S12 

Sll 

S12 

0 

0 

0 

S12 

S12 

Sll 

0 

0 

0 

0 

0 

0 

S44 

0 

0 

0 

0 

0 

0 

S44 

0 

0 

0 

0 

0 

0 

S44 

0.768 

-0.214 

-0.214 

0 

0 

0 

-0.214 

0.768 

-0.214 

0 

0 

0 

-0.214 
0 

-0.214 
0 

0.768 
0 

0 
1.256 

0 
0 

0 
0 

X  10'^^  Pa"^ 

0 

0 

0 

0 

1.256 

0 

0 

0 

0 

0 

0 

1.256 

The  modulus  of  elasticity  for  a  cubic  system  in  any  given  <  hkl  >  direction  may 
be   calculated   by   the   foUowmg   equation   in   terms    of  three   independent   elastic 
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(compliance)  constants,  Sy,  and  the  direction  cosines  of  the  crystallographic  direction 
under  study: 


1 


=  Sl] 


{^u-^n)-^H4](^p2+^p3  +^1^3) 


(C.l) 


where  /i,  I2,  l^  =  directions  cosines.  The  elastic  modulus  is  calculated  to  be  130  GPa 
in  the  <100>  direction,  169  GPa  in  the  <110>  direction,  and  189  GPa  in  the 
<  1 1 1  >  direction. 
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